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I .  INTRODUCTION 


A.  THE  MAIN  IDEA 

Image  processing  by  nonoptical  means  has  received  exten¬ 
sive  attention  in  the  last  few  years.  Several  books  and  many 
papers  have  been  published  that  have  established  nonoptical 
image  processing  as  a  viable  area  of  research.  A  large 
portion  of  this  research  emphasizes  the  linear  processing 
of  images  for  two  main  reasons:  1)  Many  image  processing 
tasks  are  linear  in  nature.  These  tasks  include  image  enhance¬ 
ment,  image  restoration,  picture  coding,  linear  pattern 
recognition,  arid  TV  bandwidth  reduction.  2)  There  are  many 
known  linear  techniques  that  may  be  brought  to  bear  in  the 
treatment  of  linear  image  processing.  These  techniques  include 
transform  theory,  matrix  theory,  filtering,  signal  modeling, 
etc.  Several  common  operations  involved  in  image  processing 
include  transfer  function  concepts,  partial  difference  (recur¬ 
sive)  equations,  and  convolution  summations.  For  example, 
Vander  Lugt  [Refs.  1,2]  has  presented  an  extensive  development 
of  linear  optics  based  on  transfer  functions.  The  transfer 
functions  relate  the  two-dimensional  Fourier  transform  of  an 
output  image  to  that  of  the  input  image.  Complex  optical 
systems  are  easily  described  by  combinations  of  transfer 
functions  that  correspond  to  individual  components  of  the 
optical  system. 


Partial  difference  equations  are  used  by  Habibi  [Ref.  3] 
to  describe  a  model  for  estimating  images  corrupted  by  noise. 
The  model  corresponds  to  a  two-dimensional  extension  of  Kalman 
filters.  Convolution  summations  are  discussed  by  Fryer  and 
Richmond  [Ref.  4]  in  work  that  involves  simplifying  a  two- 
dimensional  filter  to  a  single  dimensional  filter. 

The  time-discrete  state-space  model  offers  great 
utility  in  the  formulation  and  analysis  of  linear  sys¬ 
tems.  Linear  systems  that  are  described  by  transfer  functions, 
difference  equations,  or  convolution  summations  are  formulated 
into  a  state-space  representation.  Once  formulated,  many 
known  techniques  may  be  applied  to  systematically  analyze 
the  model.  Consequently,  the  state  space  model  is  a  general 
and  powerful  tool  that  is  used  to  unify  the  research  and  the 
study  of  time-discrete  linear  systems. 

This  thesis  develops  the  discrete  model  of  Roesser  [Ref. 

5]  for  linear  image  processing  which  closely  parallels 
the  well-known  state  space  model  for  time-discrete  systems. 
Because  it  is  parallel,  many  of  the  concepts  that  are  known 
for  the  temporal  model  may  be  carried  over  to  the  spatial 
model.  This  is  done  by  generalizing  from  a  single  coordinate 
in  time  to  two  coordinates  in  space.  The  spatial  model  will 
hopefully  have  some  of  the  same  utility  for  the  study  of  two- 
dimensional  linear  systems  as  the  temporal  model  for  one¬ 
dimensional  linear  systems  [Ref.  3].  However,  not  all  of  the 
properties  of  one-dimensional  systems  carry  over  into  the  multi 


dimensional  case. 


One  of  the  fundamental  problems  involved  with  recursive  2- 
Dimensional  systems  is  that  the  order  of  the  system  (recursive 


memory)  is  not  the  same  as  the  number  of  initial  conditions 
(boundary  conditions) .  In  one-dimensional  systems  these  are  the 
same.  Temporal  systems  are  inherently  nonanticipatory  and  are 
often  treated  as  such  for  the  sake  of  physical  realizability 
in  real  time;  whereas  spatial  systems  do  not  have  causality  which 
is  an  inherent  limitation.  That  is,  an  image  processor  may  have 
right  to  left  dependency  as  well  as  left  to  right  dependency. 
Finally  it  is  noted  that  stability  criteria  in  one-dimensional 
recursive  systems  become  much  more  difficult  when  carried  over 
to  the  multidimensional  case. 

Causality  is  built  into  the  temporal  state-space  model  if 
an  initial  state  is  assumed  to  be  fully  specified.  In  order  to 
establish  a  close  parallel  for  the  spatial  model,  the  same 
built-in  causality  will  be  intentionally  assumed,  despite  the 
fact  that  causality  is  not  necessary  for  physical  realizability 
in  real  space.  Such  an  image  processor  is  said  to  be  unilateral. 
If  the  constraint  of  causality  is  removed,  then  the  image  proces¬ 
sor  is  said  to  be  bilateral  [Ref.  5].  Concepts  that  are 
developed  in  this  thesis  for  the  latter  case  are: 

1)  Formulation  of  the  state  space  model  of  Roesser.  [Ref.  5] 

2)  The  definition  of  state  transition  matrix. 

3)  A  resulting  computer  program  based  on  the  above  model. 

4)  An  investigation  of  the  class  of  2-Dimensional  transfer 
functions  defined  by  this  model. 


5)  Derivation  of  a  general  response  formula. 


6)  Extension  of  Roesser's  model  of  state  variable  equations 
to  encompass  a  larger  class  of  transfer  functions. 

7)  Adaptation  of  the  1-D  "SSPACK"  program  to  produce  2-D  data. 

B.  STATE  SPACE  REPRESENTATION 

Toward  the  end  of  the  1950s,  the  concept  of  representing 
a  discrete  system  by  a  set  of  first-order  difference  equations 
became  a  standard  tool  of  the  research  engineer.  These  tech¬ 
niques  have  since  become  generally  known  as  state-space  repre¬ 
sentations.  Such  representations  have  become  increasingly 
important  during  the  intervening  years  because_ they  often  allow 
one  to  carry  out  a  meaningful  system  design  entirely  in  the 
discrete-time  domain  (in  comparison  to  popular  Z-transform 
methods).  That  this  is  important  follows  basically  from  these 
factors : 

1.  The  system  may  be  nonlinear  so  that  transformation 
methods  are  not  directly  applicable. 

2.  Time-domain  concepts  often  give  one  a  better  insight 
into  the  analysis  and  synthesis  of  the  system  (fre¬ 
quently  with  the  aid  of  a  digital  computer) . 

3.  Cases  in  which  the  initial  conditions  are  non-zero  may 
be  handled  straightforwardly. 

A  state  space  representation  of  a  system  differs  from  the 

conventional  representation.  In  a  conventional  representation 

only  the  relationships  between  the  input  and  output  signals  need 

be  known.  On  the  other  hand,  the  state-space  representation 

gives  a  total  description  of  both  the  internal  as  well  as  the 

external  signals  of  a  system. 
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C.  STATE-VARIABLE  REALIZATIONS — THE  CONCEPT  OF  STATE  — 

In  1-D  linear  systems  theory  and  control  theory,  the 
concept  of  a  filter  state  has  played  an  important  role. 

Basically  the  filter  state  at  any  point  in  time  contains 
all  the  information  necessary  to  compute  the  remainder  of  the 
filter  output  signal,  given  the  input  signal.  One  dimensional 
single-input,  single-output  filter  realizations  based  on  a 
state  variable  model  can  be  written  in  the  form: 

x(k+l)  =  Ax(k)  +  Bu(k)  (I. la) 

y  (k)  =  Cx (k )  +  Du (k)  (I. lb) 

This  form  relates  the  input  u(k)  and  the  output  y(k)  through 
a  state  vector  x(k).  The  state  vector  evolves  in  time  accord¬ 
ing  to  equation  (I. la).  The  matrices  A,  B,  and  C  and  lxl  matrix 
D  govern  the  exact  form  of  the  input-output  relationship. 

(In  general  these  matrices  may  vary  with  the  index  (k)  and 
the  input  and  output  signals  may  be  vectors  as  well.)  Quite 
often  the  components  of  the  state  vector  are  taken  to  be  the 
constants  of  the  Z  ^  delay  operators  in  a  flowgraph  represen¬ 
tation  of  the  1-D  filter. 

A  classic  problem  in  state-variable  theory  representation  is  to 
find  the  matrices  A,  B,Cand  D  which  will  realize  a  particular  system 
function  H(z)  with  a  minimum  number  of  state  variables.  A 
similar  approach  may  be  taken  to  develop  a  2-D  state-variable 


model . 


A  2-D  discrete  system  may  be  defined  as  a  mathematical 
abstraction  which  utilizes  three  types  of  variables  to  repre¬ 
sent  or  model  the  dynamics  of  a  discrete-time  process.  The 
three  variables  are  called  the  input,  the  output,  and  the 
state  variable.  The  input  variables  u(i,j),  serve  as  external 
forces  which  influence  the  dynamics  or  motion  of  the  system. 

The  output  variables  y(i,j)  are  the  characteristic  variables 
which  are  directly  observable  (measurable)  by  the  external 
observer.  The  state  variables  x(i,j)  characterize  the  internal 
dynamics  of  the  system  and  are  to  be  selected  according  to  the 
following  rule. 

These  variables  are  formulated  in  such  a  manner  that,  if 
one  knows  the  values  of  the  present  state  variables  x(i,j)  along 
with  the  values  of  the  input  variables  u(i,j)  then  the  output 
variables  y(i,j)  and  the  next  state  variables  x(i,j)  are  com¬ 
pletely  determined.  Moreover,  the  number  of  state  variables 
used  in  a  state-space  representation  must  be  minimized.  A 
state-space  representation  may  be  visualized  in  block  diagram 
form,  as  shown  below. 


INTERNAL 

STATES 

X,(U  )  , 


Figure  1.1 


In  Figure  1. 1 ,  m- inputs ,  p -outputs  and  n-state  variables  are  repre¬ 
sented.  However,  we  will  be  mainly  interested  in  those  systems  which 
have  one  input  (m  =  1)  and  one  output  (p  =  1) .  It  is  important 
to  note  that  the  input  and  output  variables  appear  external  to 
the  system,  while  the  state  variables  are  generally  internal. 

The  different  input  variables  will  be  represented  by  the 
input  vector  u(i,j)  where, 

u(i, j) 


ux (i,  j) 
u2 (i, j) 


the  output  vector  y(i,j)  where. 


Y(i»j)  = 


and  the  state  vector  x(i,j)  where, 


x(i,j) 


(i,  j) 
*2  (i,  j) 


For  a  given  process  the  state  space  representation  is  not  unique. 
However  all  such  representations  have  one  characteristic  in 
common  for  a  given  system,  namely  the  number  of  elements  n  is 
referred  to  as  the  order  of  the  system. 


II.  ROESSER'S  STATE-SPACE  MODEL 

A.  THE  FRAMEWORK 

An  image  is  a  generalization  of  a  temporal  signal,  in  that 
it  is  defined  over  two  spatial  dimensions  instead  of  a  single 
temporal  dimension.  Consequently,  two  space  coordinates  i  and 
j  take  the  place  of  time,  t.  Also,  two-state  sets  are  intro¬ 
duced  to  replace  the  single-state  set.  The  following  defini¬ 
tions  are  made  by  the  model: 

i  An  integer-valued  vertical  coordinate; 

j  An  integer-valued  horizontal  coordinate; 

{R}  A  set  of  n^  real  vectors  which  convey  information 
horizontally ; 

{S}  A  set  of  n2  real  vectors  which  convey  information 
vertically; 

{u}  A  set  of  m  real  vectors  that  act  as  inputs; 

{y}  A  set  of  p  real  vectors  that  act  as  outputs. 

A  specific  image  processor  is  then  defined  as  6-tuple 

< {R} , {S } , {u} , {y } , f ,g>  , 

where  f  is  the  next  state  function: 

f:  {{R},{S},{u}  +  {{R},{S}} 


and  y  is  the  output  function 


Now  since  f  and  g  are  to  be  linear  functions,  they  may  be 
represented  by  the  following  matrix  equations: 

R ( i+1 , j )  =  A1R(i,j)  +  A2S(i,j)  +  B1u(i,j) 

S ( i , j+1)  =  A^R ( i / j )  +  A^S ( i , j )  +  B2u(i,j)  (II. 1) 

y ( i , j )  =  C1R(i,j)  +  C2S(i,j)  +  Du(i,j)  i,j  _>  0 

A^ ,  A2 ,  A^,  A^,  B1  B2 ,  C^,  C2 ,  D  are  matrices  of  appropriate 

dimensions.  Boundary  conditions  R(0,j)  and  S(i,0)  and  also  the 
input  u(i,j)  are  externally  specified.  In  the  next  section  a 
computational  rule  is  obtained  that  uniquely  determines  the 
states  R(i,j)  and  S(i,j)  and  also  the  output  y(i,j)  (for  i,j  ^  0) 
from  the  boundary  conditions  (such  as  all  zero).  The  equations 
produce  a  set  of  output  vectors  from  the  input  vectors. 

This  formulation  is  general  so  that  any  discrete  linear 
image  process  may  be  so  represented.  Notation  is  condensed 
somewhat  by  introducing  the  following  matrices  and  vectors: 


T'(i,j)  =  AT ( i , j )  +  Bu ( i , j ) 

y ( i  ,  j )  =  CT ( i , j )  +  Du ( i , j ) 

B.  GENERAL  RESPONSE  FORMULA 

A  state-transition  matrix  A  is  defined  as  follows: 


Then  exponentiation  A  ,J 


is  defined  as, 


A1,j  =  A1,°A1"1,j  +  A0,1A1,j_1  (i , j )  >  (0,0) 


A0,0  =  I  ;  A-1'3  =  A1'"15  =  0  for  j  >  1,  i  >  1 


Examination  of  this  definition  bears  out  that  it  is  an  effective 
recursive  definition  of  A1 ' ^  for  integer  values  of  i  and  j such  that 
either  i  >  0  or  j  >  0  or  (i,j)  =  (0,0).  It  parallels  the  defini¬ 
tion  of  the  time-discrete  state-transition  matrix  Afc  =  A  At  \ 

It  now  remains  to  be  shown  that  this  state  transition  matrix 
A1'-1  may  be  used  in  expressions  for  the  response  of  the  model 
in  terms  of  the  inputs  and  boundary  conditions.  The  term 
boundary  conditions  is  used  here  to  refer  to  the  states  along 
the  edges  of  the  model.  Specifically,  the  set  of  boundary  condi¬ 
tions  consist  of  R ( 0 , j )  for  j  ^  0  and  S ( i , 0 )  for  i  >  0. 


C.  CHARACTERISTIC  FUNCTION  OF  A  MATRIX 

If  the  primary  inputs  and  outputs  are  dropped  in  the  model 
equations  (II. 1),  a  representation  arises  for  the  state  behavior 
of  the  system  having  the  form 


R ( i+1 , j )  =  A1R(i,j)  +  A2S(i,j) 

S ( i , j+1)  =  A^R ( i , j )  +  A^S { i , j ) 


(II.  2) 


These  equations  are  useful  in  the  development  of  a  form  for  a 
two-dimensional  characteristic  matrix  of  A.  Operators  are 
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first  introduced  that  advance  a  particular  coordinate  of 
their  operand. 

Definition:  Let  E  be  an  operator  that  has  the  effect  of  ad¬ 

vancing  the  vertical  coordinate  or  the  first  subscipt  of  the 
function  upon  which  it  is  operating.  Likewise,  let  F  be  an  opera¬ 
tor  that  has  the  effect  of  advancing  the  horizontal  coordinate  or 
second  subscript  of  the  function  upon  which  it  is  operating. 

The  effect  of  these  operators  on  the  state  vectors  is: 

R ( i+1 , j )  =  ER ( i , j ) 

S (i, j+1)  =  FS (i , j ) 

The  state  equations  can  be  rewritten  using  these  advance 
operators . 

(EI-Aj.)  R(i,  j)  -  A2S(i,j)  =  0 

-A3R(i,j)  +  (FI-A4)S(i, j)  =  0 

These  equations  are  equivalently  represented  in  the  overall 
matrix  form. 


The  above  equation  represents  a  system  of  homogeneous  equations  in 
the  elements  of  T(i,j).  If  the  system  is  to  have  a  non-trivial 
solution  for  T(i,j)  then  the  transformation  represented  by  the 
matrix  must  be  singular.  The  above  matrix  is  said  to  be  the  two- 
dimensonal  characteristic  matrix  of  the  partitioned  matrix  A,  where 


■  n  w*  v  luiiifv  u  wj  ww  w».J! 


ffiwrnmrr 


A  = 


A1  A2 
A3  A4 


i  !  o 

0,1 

0  1  0 

I 

1  o  1 

r.-f 

and  I  = 

The  characteristic  matrix  of  A  is  denoted  cm (A)  and  may  be 
represented  as 

cm  (A)  =  El1,0  +  FI0,1  -  A 

where , 


.1,0 


Now  since  cm (A)  must  be  singular,  its  determinant  must  be  equal 
to  zero.  | cm (A) |  =  0.  If  E  and  F  are  replaced  in  the  above  by 
general  indeterminates  x  and  y  respectively,  the  result  is  an 
expression  called  the  two-dimensional  characteristic  equation 
for  A.  The  determinant  of  cm (A),  and  x  and  y  replacing  E  and 
F,  is  called  the  two-dimensional  characteristic  function  of  the 
matrix  and  is  denoted  by 

|  cm  (A)  |  =  f  (x,y)  =  0 

f(x,y)  will  be  a  monic  polynomial  in  x  and  y  with  degree  n^ 
in  x,  and  degree  n2  in  y,  where  n^  is  the  dimension  of  R  and 
n2  is  the  dimension  of  S.  f(x,y)  has  the  form 


f (x,y)  = 


l  l 

(0,0)  <_(i,  j)  <_(nlfn2) 


i  ] 

a  .  .  x  y 

i,D  * 


where  a.  •  denotes  elements  of  A  and  a  =  1. 

i , 3  n^ , n2 
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D.  CIRCUIT  ELEMENTS  AND  THEIR  REALIZATION 


Let  us  consider  the  single  2-D  HR  filter  transfer  function 


given  by: 


H  ( z! ' z2 ) 


b00+b10Zl  +b01Z21+bllZllz21+b21Zl  Z2 


1-a10Zl  _a01Z2  ~allZl  Z2  "a21Zl  Z2 


B(z1,22) 

1-A(z1,z2) 


A  simple  block  diagram  for  H(z^/Z2)  follows. 


uCmj 


8  (2., 10 


y(U) 


WllJ) 


A 


3 


Figure  2.1 

The  input  signal  u(i,j)  flows  through  a  filter  corresponding 

to  the  numerator  transfer  function  B(z^,z2).  The  resulting 

signal  is  added  to  the  signal-w ( i , j )  to  produce  the  output 

signal  y(i,j).  The  denominator  transfer  function  1-A(z^,z2) 

is  realized  by  the  feedback  loop  containing  A(z^,z2). 

Since  we  are  dealing  with  two  dimensions,  there  are  two 

fundamental  shift  operators  which  may  occur  along  a  signal  flow 

path,  the  horizontal  shift  operator  indicated  by  z^1  and  the 

vertical  shift  indicated  by  z3  [we  shall  omit  from  consideration 

>  z. 

the  inverse  shift  operators  z^  and  z2] .  In  most  cases  of  prac¬ 
tical  interest  they  can  be  eliminated  by  multiplying  both  the 
numerator  and  denominator  polynomials  of  H(z^,z2)  by  the  appro¬ 
priate  powers  of  and  z^ .  Let  us  look  at  a  signal  flowgraph 

representing  the  numerator  polynomial: 


Figure  2.2 


Note  the  chain  of  two  z^  operators  descending  on  the  left  and 
the  single  z^~  operator  ascending  on  the  right.  The  nodes 
along  these  two  vertical  paths  are  connected  by  branches  with 
the  appropriate  gains.  If  we  label  the  nodes  in  both  z~l 
chains  and  the  z^~  chain  0,1,2  and  so  on,  from  the  top  down, 
the  ith  node  in  the  z^1  chain  is  connected  to  the  jth  node 
in  the  z^-  chain  by  a  branch  with  a  gain  factor  of  b^ ^ . 

Similarly  the  signal  flowgraph  for  the  polynomial  A{z^,z^) 
is  shown  in  Figure  2.3. 


OWTfWT 


Figure  2.3 


Since  there  is  no  a qq  term,  there  is  no  direct  connection  be¬ 
tween  the  input  and  output  nodes  of  this  signal  flowgraph.  Thus 
any  path  from  the  input  node  to  the  output  node  will  encounter 
at  least  one  z^  or  shift  operator. 

At  this  point  it  is  appropriate  to  discuss  realizations  for 
the  two  shift  operators  z^1  and  ^2  '  At  their  simplest  level, 
the  shift  operators  merely  select  the  "previous"  S-tuple  value 
in  the  horizontal  or  vertical  direction.  When  the  input  to  a 
z11  operator  is  the  S-tuple  u(i,j)  the  output  will  be  R(i-l,j). 
Similarly  for  a  operator  the  output  will  be  S  ( i , j  —  1 )  when 

the  input  is  R(i,j)  or  S(i,j).  Consequently  a  realization  of 
either  shift  operator  must  embody  the  appropriate  amount  of  memory 
to  retain  the  "previous"  S-tuple  in  the  appropriate  direction. 

Interestingly  enough,  in  the  more  general  case  where  the  numer 
ator  and  denominator  polynomials  are  considered  jointly,  the  state 
variable  realizations  based  on  conventional  signal  flowgraphs  may 
not  be  minimal  in  the  sense  that  the  transfer  furnction  can  be 
realized  with  fewer  coefficients.  Consider, 


H  ( zi ' Z2 ) 


b10  1  b01  2  bll  1  Z2 

:  =t - =r - -T"-r 

1_a10Zl  _a01Z2  -allZl  Z2 


(II. 3) 


The  corresponding  signal  flow  representation  is  shown  in  Figure 


2.4  below: 


E.  ANALYSIS  OF  ROESSER'S  MODEL 

Recalling  from  page  14  the  equations  of  the  model  are 


R  (  i+1 r j )  =  A1R(i,j)  +  A2S(i,j)  +  B1u(i,j) 
S ( i / j+1)  =  A^R ( i / j )  +  A4S(i,j)  +  B2u(i,j) 
y  ( i / j )  =  C1R(i/j)  +  C2S(i,j)  +  Du ( i , j ) 


fc: 


■ty'v; 


p 


A A2 ,  A^ ,  A^ ,  B^ ,  B2/  C1 ,  C2 ,  D  are  scalars  or  matrices  of 
appropriate  dimensions. 


R ( i+1 , j  ) 
S  (i ,  j  +  1) 


A1  A2 
A3  a4 


R(i,  j) 
S  ( i ,  j  ) 


+ 


B. 


B, 


u(i, j) 


y(i,j)  =  tc1  c2] 


R(i, j) 

S(i,  j) 


+  Du ( i  ,  j  ) 


R( i+1 ,  j  )  =  A}R(i,j)  +  A2S(i,j)  +  B1u(i/j) 
S ( i+1 , j )  =  A^R { i ,  j  )  +  A^S ( i , j )  +  B2u(i,j) 


(II. 4) 


(II. 5) 


And  taking  Z  transforms: 


z1R(z1,z2)  =  A1R(z1,z2)  +  Z2S(z1,z2)  +  b1u(z1,z2) 


=  A^R  (  z]_ '  Z2  )  +  A^S  (  ,  z2  )  +  B2u(z1,z2) 


22S(z1'V 


22S(z-l,z2)  -  A3R(z1,z2)  -  A4S(z1,z2)  =  B2u(z 

Dr 

r(z1/Z2)  [zx  -A11  -  A2S(z1,z2)  =  B1u(z1,z 

R(z1/Z2) [— A3 ]  -  [z2-A] S (z1/Z2)  =  B2u(z1,z 


+  Du  (z,  ,  z_ ) 


0 


(ii. 7; 


The  submatrix  Z-^  is  simply  z^  times  an  identity  matrix  of  the 

appropriate  size.  Similarly  z2  is  z2  times  an, identity  matrix. 

The  objective  of  the  state  variable  realization  procedure  is 

to  find  the  matrices  A,  B,  C,  and  D  which  yields  an  F(z  ,z  ) 

-L  2.  • 

that  equals  or  approximates  a  desired  system  function  H(z^,z2). 

In  essence,  the  equations  of  Roesser  represent  an  implementation 
for  which  a  design  algorithm  must  be  found.  One  choice  for 
the  state  variables  is  the  output  signals  from  the  shift 
operators . 

Thus  R(i,j)  is  a  vector  containing  the  output  signals 
from  the  z^1  operators  and  S(i,j)  contains  the  output  signals 
from  the  z2^  operators.  (Note  that  the  output  signal  of  a 
shift  operator  signal  path  is  not  necessarily  the  same  as  the 
nodal  signal  at  the  node  to  which  the  signal  path  points.) 

If  a  state  variable  corresponds  to  the  output  of  a  shift  operator, 
the  next  value  of  that  state  variable  must  correspond  to  the 
input  of  the  shift  operator.  To  obtain  the  submatrices  , 

A9 ,  A^ ,  A^  in  equations  of  Roesser,  we  write  the  input  signal 
of  each  shift  operator  in  terms  of  the  outputs  of  all  the 


shift  operators,  taking  care  to  include  all  shift-free  paths 
from  ouput  to  input  (see  the  following  flowgraph) . 


CC1  C23 


[Cl  c2] 


(22-A4)Bl 

(Z2'A1)  (A2-Z4) -A2A3 

A3B1 _ 

( Z1-A1) (Z2-A4) -A2A3 


(Z2'A4)B1  +  A2B2 

(Z2-A1) (Z2-A4) -A2A3 
A3B1  +  (zi_a1)B2 

(zi-ai) (z2-a4) -a2a3 


A2B2 _ 

(Z1~A1) (Z2-A4) -A2A3 

(2rAi)B2 

(Z1-A1) (Z2-A4) -A2A3 


H  ^  z]_ '  z2 ) 


Ci  (  z2~A4  )  Bi  +  C-LA2B2  +  C2A.3B1  +  C2  (  Z 1“ Ai )  B 2 


(Zi-A1) (Z2-A4) -A2A3 


h(z1'Z2) 


ClBlZ2~ClBlA4+ClA2B2+C2A3Bl+C2B22l"’C2B2Al 
^1Z2  _A421  ~  A1Z2  _  A2A3  +  A1A4 


(ClA2B2+C2A3Bl~C2B2Al"ClBlA4) + (C2B2Z1+C1B1Z2) 


( A1A4~A2A3  ~ A4Z1^A122^212 2 


Equating  equation  (II. 8)  with  (II. 3)  on  page  21  yields 


(II. 8) 


101  012 

a 

n~i  2 

For  this  example,  =  z^,  Z^  = 

Z2 ' 

all  of  the 

coefficients  on 

the  left  hand  side  are  scalars. 

Equation  terms 

of  equal  powers 

of  z1  and  z 2, 

C-,  A_B_+C-,A_B..  -C0B-A-C.B_A. 
12  2  231  221  114 

= 

bu  -  0  ' 

(II. 9) 

C2B2. 

= 

bio 

(II. 10) 

C1B1 

= 

boi 

(II. 11) 

A1A4  ~  A2A3 

= 

1 

(11.12) 

A4 

= 

aio 

(11.13) 

A1 

= 

a0l 

(11.14) 

all 

= 

-1 

(11.15) 

From  these  equations,  assuming  that 

B.  =  B,  =  1 

,  it  follows  that: 

From  Equation  (11-12) : 


A1A4  " 


A2A3 


=  1  =  -a 


11 


-A2A3  =  'aU  -  a1A4 
A2A3  =  all  +  a10a01 

Let  A2  and  take  on  particular  values  p  and  q  respectively, 

A3  =  q  As  "  P 

or 


pq  =  a,,  +  aina 


11  “10“01 


From  Equation  (II-6)  : 


C.  AnB_  +  C-,A_B.  —  C-BnA,  -  C.B-A  —  b.  . 
122  231  221  114  11 


or. 


b01P  +  bioq  “  bioaoi  "  boiaio  ~  bll  0 


Substituting  Equation  (11.16)  into  Equation  (11.17): 


(II 


(II 


an+aioaoi 

b01  q  +  b10q  “  b10a0l  "  b0lai0  "  bn  “  0 


.16) 


•  17) 


0  •  (11.18) 


orii  “ono  or 


The  results  are  just  the  same  as  in  [Ref.  8] .  After  the  com¬ 
parison  between  Roesser's  model  and  the  2-D  IIR  filter, 
described  by  Equation  (II. 3),  we  have: 


p  ;  pq  =  an  =  a10aQ1 


^ “  ^10aQl+^0iai0+^ll 1  3-1  1  +^n  1  3-1  1  ^ 


oi  ii  orio  or 


(11.19) 


The  foregoing  equations  relate  the  coefficients  of  the  2-D 
transfer  function  to  the  terms  of  the  system  matrices  of  the 
Roesser  model,  Equation  (II. 1). 

Rung  et  al.  [Ref.  2]  have  shown  that  the  following  state 
variable  equations,  which  use  only  two  shift  operators,  will 
also  realize  .  For  the  foregoing  example, 


mmmm 
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We  can  construct  a  signal  flowgraph  with  only  two  shift 
operators.  It  is  an  equivalent  figure  to  that  on  page  25. 

Rung  et  al .  [Ref.  2]  have  also  shown  that  state-variable 
realizations  of  the  form  of  the  equations  above  may  be 
generalized  for  any  system  function  H which  satisfies 
the  following  three  conditions: 

°<.o 


1)  The  constant  term  in  the  numerator,  bQQ  =  0,  must  be 
zero. 

2)  The  largest  powers  of  z-j_\  in  the  numerator  and 
denominator  polynomials,  must  be  equal,  and 

3)  The  largest  powers  of  in  the  numerator  and 

denominator  polynomials  must  be  equal. 

There  is  one  potential  difficulty  with  state  variable  reali¬ 
zations  of  this  type.  The  nonlinear  equations  defining  p  and 
q  may  result  in  complex  values  for  these  constants.  For 
example,  when  b^Q  =  bQ1  =  1,  b^  =  0,  a1Q  =  aQ1  =  2  and 
a...  =1,  we  get  p  =  q*  =  2  ±  j. 


-.'V.  s’  •.VV/.V.W 


III.  THE  PROGRAM  OF  ROESSER’S  EQUATIONS 

WITH  SCALAR  COEFFICIENTS  (FIRST  ORDER) 


A.  AN  EXAMPLE 

For  a  4  x 4  data  field  the  S  and  R  matrices  are  indexed 
as  follows: 


j  j 


i 

'i 


S  matrix 

For  4x4  Matrices 

The  Initial  Conditions  are  given,  by  the  values 
R(l,l) ,  R(2 , 1) ,  R  ( 3 , 1) ,  R ( 4 , 1) 

S  ( 1 , 1 )  ,  S  ( 1 , 2  )  ,  S  ( 1 , 3  )  ,  S  ( 1 , 4  ) 

The  2-D  state  variable  equations  can  be  written  as: 


1.1  1,2  1,3  1,4 

2.1  2,2  2,3  2,4 

3.1  3,2  3,3  3,4 

4.1  4,2  4,3  4,4 

R  matrix 


R  (  i  +  1 , j ) 

S  ( i ,  j+1) 

y(i,j)  = 


=  A-,  R  ( i  ,  j  )  +  A-S ( i , j )  +  B. u ( i  ,  j  ) 


=  A^R ( i , j )  +  A4S(i,j)  +  B2u(i,j 


tct  c2] 


R  (  i  ,  j  ) 
S (i, j) 


The  input  2-D  field  is  taken  to  be, 


u  (  i , j )  =  1 ,  for  i  =  j  =  1 


0 


otherwise . 


The  output  data  field  is  indexed  as 


_ __J _ 

1.1  1,2  1,3  1,4 

2.1  2,2  2,3  2,4 

i 

,  3,1  3,2  3,3  3,3 

4.1  4,2  4,3  4,4 

Y  output  matrix 

B.  THE  2-D  FOURIER  TRANSFORM 

The  2-D  discrete  Fourier  transform  Y(m,n)  q f  the  output 
y(i,j)  can  be  written  as, 


M-l  N-l  -j2TT—  -j2TT— 

Y  (m,n)  =  l  l  y  ( £ ,  k)  e  M  e  N 

£=0  k=0 


or  for  convenience. 


Y  (m,n) 


M  N  j  2n 

l  j  y  (  £  ,  k)  e 
£=1  k=l 


(£-1) (m-l ; 
M 


—  j  2  TT 


(k-1) (n-l) 
N 


Y(m,n):  2-D  D.F.T.  (y(i,j) } 

M  *N:  The  dimension  of  the  given  data  y(£,k)  and  D.F.T. 
Y(m,n)  also. 

y(£,k):  Given  data  (The  output  as  described  above). 

To  develop  the  D.F.T.  for  two-dimensional  signals  we 
consider  a  finite  area  sequence  y(£,k)  which  is  zero  outside 
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*.< 


.a 


it  is  of  area 


R 


w. 


£ 


fS 


‘i 


is 


K 


the  interval  0  <_  £  £  M-l,  0  £  k  £  N-l,  i.e., 
(M,N)  and  construct  the  periodic  sequence: 


y(£,k)  =  y  [((£))  M  (( k ))  N )  ] 


The  original  sequence  y(£,k)  is  recovered  by  extracting  one 
period  of  y(£,k),  i.e.. 


y(£,k)  =  y  (£,k)  RM  N(«.,k) 


RM,N(£,k) 


1,  0  <  £  £  M-l ,  0  <„  k  £  N-l 

0,  otherwise  . 


We  then  define  the  discrete  Fourier  transform  of  y(i,k)  to 
correspond  to  the  Fourier  series  coefficients  of  y(£,k). 
However,  just  as  we  did  with  one-dimensional  sequences,  we 
will  maintain  the  duality  between  the  time  and  frequency 
domains  by  interpreting  the  D.F.T.  coefficients  to  also  be  a 
finite  2-D  sequence.  Thus  with  Y(m,n)  denoting  the  D.F.T. 
of  y(£,k),  we  can  write 


M- 1  N- 1 


Y (m,n) 


. ~  £m  . „  kn 
-1)2^—  -  D  2tt— — 

l  l  yU,k)e  e  R  «(m,n) 

£=0  k=0 


or 


y  U,k) 


M-l  N-l 

=  i  I  I  Y(m,n)  e 


j  2tt 


£m 


m=0  y=0 


•  p  kn 

M  eD  N  R  (?,k) 
M ,  N 
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or, 

M  N  -j2TT— ----  -j27T-— -1)  (n~1) 

Y(m,n)  =  l  l  y.  .(£,k)e  N  e  N 

£=1  k=l  1,11 

As  an  example,  consider  the  case  forNM  =  N  =  5. 

Given  2-D  Data  Sequence 

1.1  1,2  1,3  1,4  1,5 

2.1  2,2  2,3  2,4  2,5 

y(£,k)  =  3,1  3,2  3,3  3,4  3,5 

i  4,1  4,2  4,3  4,4  4,5 

5.1  5,2  5,3  5,4  5,5 

j 

Then , 

y (1,1)  +  y(l, 2)  +  y (1,3)  +  y(l,4)  +  y (1,5) 

+y(2,l)  +  y (2,2)  +  y(2,3)  +  y(2,4)  +  y(2,5) 

Y(l,l)  =  +y(3,l)  +  y (3,2)  +  y(3,3)  +  y(3,4)  +  y(3,5) 

m=1,n=1  +y(4,l)  +  y (4,2)  +  y (4,3)  +  y(4,4)  +  y(4,5) 

+y (5,1)  +  y (5,2)  +  y(5,3)  +  y(5,4)  +  y(5,5) 


Matrix  5x5 

M=5  N=5 

£=1,2, 3, 4, 5  k=l , 2 , 3 , 4 , 5 


4(-jc)  8(-j£)  12  (~jf  16  (- j- 

+y(5,l)  +  y (5,2)e  3  +  y(5,3)e  3  +  y(5,4)e  3  +  y(5,5)e 


In  Appendix  A  we  give  a  listing  of  the  programs  that 
been  written  to  generate  y(i,j)  and  Y(m,n). 

C.  NUMERICAL  EXAMPLES 

Three  numerical  examples  which  depend  on  Equation  (I] 
are  used  to  demonstrate  the  program  in  Appendix  A. 

First  example: 


yields 


h(z1#z2: 


•  5  (z11  +  z”1) 

:  r^r 

1  -  .  2Zi  -  .  3z^ 


lll 


Bi 


B, 


D 


=  a 

=  a 

-  b 

=  b 

=  1 
=  1 

=  0 


01 


10 


10 


01 


0 

0.3 

0.2 

0.5 

0.5 


have 


.19) 


After  substitution  of  these  values  in  Eqs.  (11.16)  and  (11.18) 
we  identify 


0.5 


(III. la) 


&5I 


Second  example: 

Proceeding  in  a  similar  way  with 


H(z1,z2) 


0.25z11+0.3z21+0.'2z11z21 

1-0.125z“1-0.2z21-0.1z“1z21 


yields 


11 

=  0.1 

boo 

=  0  . 

01 

=  0.2 

bll 

=  0.2 

= 

0.125 

= 

0.83 

= 

1  or 

= 

0.15 

1  A 

=  0.125 

(III. lb) 


Third  example: 


Proceeding  in  a  similar  way  with 


0.25z~1+0.15z~1+0.72z~1z~1 

K(zl'z2}  =  - “Tf - -  -f  -] 

l-0.135z1±-0.25z  -0.15z  z 


yields 


=  0  b. 


=  0.25 


0.15 


=  0.25 


=  p  =  0.1312 


=  q  =  1.4 


(III.lc) 


=  0.135 


=  0.15 


=  0.25 


Zero  initial  conditions  were  assumed  for  all  examples.  The 


simulation  results  are  presented  in  Figures  3-1,  3-2  and  3-3 


Contour  Map  for  Figure  3-la 


Figure  3-3a.  2-D  D.F.T.  Sequence,  Y( 


Figure  3-3b.  Contour 


In  order  to  verify  the  correctness  of  the  output  produced 
by  Roesser,  the  2-D  D.F.T.  Y(m,n)  plots  for  these  examples 


were  compared  with  the  corresponding  |h(z^,z2)|.  2-D  transfer 

function  plots  |h(z^,z2) |  for  Examples  1,  2  and  3  are  shown 
in  Figs.  3-4a,b,  3-5a,b  and  3-6a,b  respectively.  The  listing 
of  a  program  used  to  generate  these  plots  can  be  found  in 
Appendix  B. 
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Transfer  Funct 
for  Example  2 


IV.  EXTENSION  OF  ROESSER'S  MODEL  TO  SECOND  AND  HIGHER  ORDERS 


A.  MINIMIZING  THE  NUMBER  OF  SHIFT  OPERATORS 

In  order  to  minimize  the  number  of  shift  operators  we  follow 
the  procedure  given  in  Kung  [Ref.  8] .  Let  us  consider  the 
simple  2-D  HR  filter  transfer  function  given  by 


H  (z! ' z2 )  = 


b00  b10  1  b01Z2  bll~l  z2  *b21Zl  22 

"  -1  -1  -1  -1  -2  -2  -1 
1-a10  1  _a0l  z  "allZl  Z2  ‘a102l  "a21Zl  Z2 


B(Zl,z2) 

1-A ( z^ , z2) 


(IV. 1) 


Our  problem  will  be  drawing  a  detailed  signal  flowgraph  for 
the  system  function  H(z^,z2).  We  can  do  this  simply  enough  by 
combining  the  flowgraphs  on  Figures  2-2  and  2-3  to  get  the 


Figure  4-1 


This  flowgraph  can  be  made  even  simpler  because  the  shift  opera¬ 
tion  is  distributive  over  addition.  We  can  combine  the  two 
operators  into  a  single  one,  yielding  the  following  flowgraph. 


Figure  4"-3 


Then,  when  we  substitute  Figures  2-2  and  2-3  for  the  blocks  as 
before,  the  two  z^  chains  will  contain  the  same  data  and  can 
be  merged  to  yield  the  signal  flowgraph  in  Figure  4-4.  This 
glowgraph  has  a  total  of  four  shift  operators,  and  it  mini¬ 
mizes  the  number  of  z^-  operators. 


N*WT 
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Figure  4-4 


Another  signal  flowgraph  that  minimizes  the  number  of  z 
operators  may  be  obtained  from  Figure  4-4  by  the  2-D  transposition 
theorem  to  obtain  a  transposed  network.  Like  i-ts  1-D  counter¬ 
part  [Ref.  9],  the  2-D  transposition  theorem  states  that  the 
transposed  network,  which  is  obtained  by  reversing  the 
directions  of  all  the  arrows  in  a  signal  flowgraph,  will  have 
the  same  system  function  as  the  original  network.  If  we 
reverse  the  direction  of  all  the  arrows  in  Figure  4-4  and  then 
redraw  the  flow  graph  with  the  input  port  on  the  left  and  the 
output  port  on  the  right,  we  get  the  flowgraph  shown  in  Figure  4-5 


OmT 


lc3 


Figure  4-5 


This  transposed  flowgraph  may  be  preferred  in  implementations 
with  limited  wordlengths  since  the  attenuation  due  to  the 
"zeros"  of  HCz^/Z^)  occurs  before  the  gain  due  to  the 
"poles"  thus  lessening  somewhat  the  possibility  of  arithmetic 
overflow  in  the  intermediate  computations. 

Using  the  notion  of  transposition  at  both  the  flowgraph 
level  and  the  block  diagram  level  (note  that  Figure  2-2  is  the 
transpose  of  Figure  2-1)  the  flowgraph  can  be  manipulated  to 
yield  a  realization  that  minimizes  the  total  number  of  shift 
operators . 

As  we  saw  earlier,  however,  a  operator  will  require 

substantially  more  storage  than  a  z^l  operator  for  a  row-by¬ 
row  ordering  of  input  samples.  Consequently,  it  may  be  more 
economical  to  minimize  not  the  total  number  of  shift  operators 
(as  in  the  1-D  case)  but  the  number  of  operators. 

If  the  filter  is  realized  by  using  a  separate  micro¬ 
processor  to  compute  samples  of  each  node  signal,  storage 
may  be  less  of  an  issue. 

In  this  case,  we  may  want  to  minimize  the  total  number  of 
nodes  in  a  flowgraph  in  order  to  reduce  the  number  of  micro¬ 
processors  in  an  implementation. 

As  digital  technology  progresses,  the  relative  costs  of 
storage,  computation,  and  interconnectivity  keep  changing.  In 
the  future  digital  systems  designers  may  have  radically  differ¬ 
ent  criteria  for  optimizing  a  filter  realization. 


I 


B.  A  SECOND  ORDER  MODEL 

Looking  at  the  flowgraph  in  Figure  4-5  and  developing  a 


state  variable  implementation  from  it,  we  shall  call  the  output 

_  -C  _ —  1  _  _  — _ -  T-v  t  Z  '  \  1  _  _ 1 _ _  i  _  —  1  1  _  _  -1 _ 


of  the  top  z^  operator  R^(i,j),  the  output  of  the  lower 
operator  R2(i,j),  the  output  of  the  left  z2^  operator 
S^(i,j)  and  the  output  of  the  right  Z21  operator  S2(i»j)  as 
indicated : 


ucu 


Figure  4-6 


H  ( z! ' z2 ■ 


b00+b10Zl1+b01Z21+bllZl  z21+b21Zl  z2 


1“a10Zl  ~a01Zl1"allZl  Z21_a20Zl  "a21Zl  Z2 


R1(i+l, j) 

R2 (i+1, j) 

S1  (i , j+1) 

S2  (i,  j  +  1) 

1  bn+boiaio  au+aoiaio 


Rx (i»  j ) 


0  b21+b0ia!0  a21+a0ia20 


|  r2 (i  '  D  ) 

S1(i,  j) 


s2  ( i » j  ) 


bio+booaio 


b00a20 


u(i, j) 


(IV. 2) 


•2535 


Y(i,j)  = 


11  0  b0l  aoi] 


RL ( i / j ) 
r2  j  ) 

S1(i, j) 
S2  (i,  j) 


+  [bnn]u(i,j)  (IV. 


Defining 


bll  =  bll  +  aioboi 


all  all  +  a10a0l 


b21  b21  +  a20b01 


a21  a21  +  a20aGl 


In  general  the  foregoing  equations  can  be  written  as 


bij  “  bii+ai0b0j 


a •  ■  =  a  .  . +a . _a„  , 

13  i]  1O  O3 


Now  we  can  give  an  expanded  version  of  (IV-2) 
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Equations  (IV. 2)  and  (IV. 3)  represent  an  algorithm  for  comput¬ 
ing  the  samples  of  the  output  signal  from  the  samples  of  the 
input  signal.  Just  as  in  the  preceding  subsection,  the  amount 
of  memory  required  to  store  the  state  variables  depends  on  the 
order  in  which  the  output  samples  are  to  be  computed.  It  is 
possible  to  envision  a  multiprocessor  architecture  for  comput¬ 
ing  equation  (IV. 4)  by  assigning  each  processor  the  responsi¬ 
bility  of  computing  the  next  value  of  a  particular  state 
variable  given  the  current  input  value  and  the  current  state- 
variable  values.  Equation  (IV. 3)  could  be  implemented  by  a 
filter  microprocessor  to  generate  the  desired  output  signal 
values . 

In  such  an  architecture,  minimization  of  the  number  of 
microprocessors  corresponds  to  the  minimization  of  the  number 
of  state  variables,  a  problem  studied  thoroughly  in  the 
literature.  Other  state-variable  forms  with  the  same  number 
of  state  variables  can  also  be  found  that  will  realize  the 
same  system  function  Htz-^fZ^)  and  may  exhibit  lower  coefficients 
of  sensitivity  or  round-off  noise  [Refs.  2,10]. 

For  the  special  case  of  "all-pole"  2-D  HR  filters,  that  is, 
filters  with  a  system  function  of  the  form: 


where  bQQ  is  a  constant  and  A(z^,z2)  is  a  2-D  polynomial,  it 
can  be  shown  that  state  variable  realizations  based  on  signal 


flowgraphs,  using  the  output  of  the  shift  operators  as  the 
state  variables,  require  the  minimum  number  of  state  variables. 
They  are  minimal  realizations  [Ref.  2] . 

From  the  above  equations  corresponding  to  the  second  order 
Roesser  model,  the  program  in  Appendix  C,  was  written.  This 
program  uses  the  values  of  coefficients  of  H(z^,z2)  as  inputs 
and  it  generates  an  output,  y(i,j).  Next,  the  program  finds 
the  2-D  Fourier  transform  of  this  output  matrix,  and  compares 
it  to  the  transfer  function 
Numerical  Example 

In  the  following  three  examples  (first  and  second  orders), 
we  use  the  coefficients  of  first  and  second  order  transfer 
functions.  We  consider  the  special  case  of  "all-pole"  2-D  HR 
filters,  i.e.,  filters  with  a  transfer  function  of  the  form: 

H ( z  z  )  —  b0Q  _  1 

“  Z1 '  Z2  "  A(Z;l,z2) 

where  bQQ  is  constant  (unity  in  our  case)  and  A(z^,z2)  is  a 
2-D  polynomial.  It  can  be  shown  that  state  variable  realiza¬ 
tions  based  on  signal  flowgraphs,  using  the  output  of  the  shift 
operators  as  the  state  variables ,  require  the  minimum  number  of 
state  variables.  They  are  minimal  realizations  [Ref.  2]. 

For  the  third  program  we  have  a  graph  for  the  case  of  a 
BP  filter. 


Example  #4 


The  |Y(m,n) |  for  this  example  is  plotted  in  Fig.  4-la.  The 
corresponding  contour  map  of  2-D  surface  is  shown  in  Figure 
4-lb. 

Example  #5 


H(z  ,2  )  =  - — - T7 - rr-r I - T2-TX  (IV. 7) 

l-0.t5z11-0. 345z21-0.125z1  z2  -0.1z1  z2 

The  2-D  D.F.T. |y(m,n)|  of  the  output  of  this  filter  is  shown 
in  Fig.  4-2a.  The  corresponding  contour  map  is  shown  in 
Fig.  4-2b. 

Example  #6 


H ( z! ' z2 ) 


-0. 125+0. 25z~1+0. 125z21-0 .125z11z21+0.125z^2z21 


1  + 


(IV. 8) 


|Y(m,n)|  for  this  examle  and  the  corresponding  contour  map  are 
shown  in  Figs.  4-3a  and  4-3b,  respectively. 

For  reasons  of  verification,  as  before,  Y(m,n)  was  com¬ 
pared  to  the  actual  transfer  function  H(u)^,u>2)  for  examples 
4,5,6.  These  transfer  functions  plots  and  the  corresponding 
contour  maps  are  shown  in  Fig.  4-4a,b,  Fig.  4-5a,b  and  Fig. 
4-6a,b  for  the  examples  4,5,6,  respectively. 

C.  EXTENSION  OF  THE  2-D  STATE  SPACE  MODELS  TO  HIGHER  ORDER 
TRANSFER  FUNCTIONS 

1 .  Introduction 

During  recent  years,  several  authors  (Attasi  [Ref.  11, 


Fozmasimi  and  Mazchesini  [Ref.  13],  Givone  and  Roesser  [Ref. 
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Figure  4-2a.  2-D  D.F.T.  Sequences  Y(m,n)  for  Example 
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Figure  4-3b.  Contour  Map  for  Figure  4-3a 


Figure  4-4a.  Transfer  Function  jn(z 
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Figure  4-5b.  Contour 
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Figure  4-6b.  Contour  Map  for  Figure  4-6a 
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17])  have  proposed  different  state  space  models  for  2-D 
systems.  They  have  also  suggested  some  extensions  of  the 
usual  1-D  notions  of  controllability,  observability,  and 
minimality  to  the  2-D  case. 

However,  these  results  are  not  quite  satisfactory. 

They  either  lack  motivation  for  the  state-space  models  intro¬ 
duced  or  the  notion  of  state-space  is  improperly  defined.  In 
Chapter  II  we  started  with  a  comparison  of  all  the  current 
models  based  on  a  practical  (circuit-oriented)  point  of 
view  and  on  a  proper  definition  of  state.  It  is  shown  that 
the  model  of  Roesser  is  the  most  satisfactory,  in  that  it  is 
also  the  most  general  since  the  Attasi  and  Fozmasimi  Mazchesimi 
models  can  be  imbedded  in  the  Givone  and  Roesser  model. 

In  Chapter  II  we  pointed  out  that  a  major  difference 
between  1-D  and  2-D  systems  is  that  in  the  2-D  case  a  global 
state  (which  preserves  all  past  information)  and  a  local  state 
(which  gives  us  the  size  of  the  recursions  of  the  2-D  filter) 
can  be  introduced. 

2 .  Extension  for  2-D  Systems 

In  [Ref.  14],  Fozmasimi  and  Mazchesimi  use  the  algebraic 
point  of  view  of  "Nerode"  equivalence.  In  this  framework,  the 
state  space  arises  from  the  factorization  of  the  2-D  input/ 
output  map.  Fozmasimi  and  Mazchesini  were  the  first  to  realize 
that  a  major  difference  between  1-D  and  2-D  systems  is  that  we 
can  introduce  a  global  state  and  a  local  state  in  the  2-D  case. 

The  global  state  (which  is  of  infinite  dimensions,  in 
general)  preserves  all  the  past  information  while  the  local 


state  gives  us  the  size  of  the  recursions  to  be  performed  at 
each  step  by  the  2-D  filter.  However,  Fozmasimi  and  Mazchesini 
failed  to  exploit  fully  the  structure  of  the  global  state  and 
its  relation  to  the  local  state,  so  that  the  state  space  model 
they  introduced  is  unsatisfactory  in  the  sense  that  what  they 
introduce  as  the  state  is  really  only  a  "partial  state"  (as 
defined  by  Wololich  [Ref.  15]  for  1-D  systems).  Indeed,  this 
partial  state  does  not  obey  a  first-order  difference  equation 
(the  notion  of  first  order  difference  equation  for  linear  sys¬ 
tems  or  partially  ordered  sets  has  been  defined  by  Mullans  and 
Elliot  in  [Ref.  1 6 ] X  Attasi's  model  suffers  from  the  same 
drawback  as  the  Fozmasimi  and  Mazchesini  one. 

On  the  other  hand,  Givone  and  Roesser  [Refs.  17,18,1] 
have  used  a  "circuit  approach"  to  the  problem  of  state  space 
realization  for  2-D  systems.  They  present  a  model  in  which 
the  local  state  is  divided  into  a  horizontal  and  a  vertical 
state  which  are  propagated,  respectively,  horizontally  and 
vertically  by  first-order  difference  equations.  From  this 
point  of  view  the  global  state  appears  as  the  boundary  condi¬ 
tion  necessary  to  propagate  the  state-space  equations. 

However,  Roesser  did  not  provide  much  motivation  for 
the  introduction  of  such  a  model  and  seemed  unaware  of  the 
full  circuit  interpretation  of  their  model  since  they  were  not 
able  to  implement  an  arbitrary  2-D  transfer  function,  say 


Mitra  et  al  gave  an  answer  in  [Ref.  19]  by  presenting 


an  implementation  method  for  2-D  transfer  functions  using  delay 
elements  z^3-  and  .  We  shall  see  below  that  this  approach  is 
consistent  with  Roesser's  model.  It  is  shown  in  [Ref.  8]  that 
Roesser's  model  appears  naturally  as  a  way  to  describe  the  local 
state  properties.  For  a  (n,m)  2-D  transfer  function, 


H(zL,z  ) 


b(Zl,z2) 
a  ( 2i ' Z2 ) 


n  m 


l  l  hijzilz23 

i=0  j=0  13  1  z 


n  m 

l  l 

i=0  j=0 


j  l  a  .  .  z  z 
,  i  ~  ~  i]l  2 


(IV. 9) 


exhibits  some  canonical  state-space  forms  (controllability, 
observability),  which  can  also  be  written  as, 


(IV. 10) 


Without  loss  of  generality,  we  can  assume  =  1  and  we 

denote 


a0 (z^  )  1  +  aQ (^  ) 

Thus,  using  1-D  realization  technique,  H(z^,z2)  of  Eq .  (IV. NO) 
can  be  used  as  shown  below  in  Fig.  4-7. 


Figure  4-7 


The  realization  is  almost  achieved:  in  addition  to  the 
n-horizontal  delay  elements,  we  need  only  m  vertical  delay 
elements  to  implement  the  feedback  gains  {a^(z  ■*"),  i  =  0,1,..., m 
and  m  other  vertical  delay  elements  to  implement  the  readout 
gains  {b^lz^),  i  =  0,1,..., m}.  Thus  the  complete  realization 
shown  in  Fig.  4-8  requires  only  n+2m  dynamic  elements.  This 
realization  is  a  standard  (canonical)  one;  its  structure  is  very 
simple  and  it  involves  only  real  gains.  Note  also  that  we 
need  fewer  dynamic  elements  than  was  suggested  by  the  imple¬ 
mentations  of  [Ref.  19]. 

As  mentioned  in  Section  (b),  circuit  implementations 
with  delay  elements  z^~  and  z^~  are  in  a  one-to-one  corres¬ 
pondence  with  state-space  models  of  Roesser's  type.  The 
outputs  of  the  z^l  delays  are  the  horizontal  states  and  the 
outputs  of  the  z^~  delays  are  the  vertical  states.  Thus  the 
implementation  of  the  following  figure  can  be  transformed 
readily  into  the  following  state-space  model. 


'  f.  C 

“si  :k:i 


i  S32  Ia  2  a22  hV2 


a02  '°l*  l°:2  [°n 


I  I  I  ' 

laOJ  {^7t  la3. 


2.0,0 


K\C.'0') 


s.coo 


_  3:a  ■  °-°  y^,  (jjj 

Figure  4-8 


R  (  i  +  1 , j ) 

R(i,  j ) 

S1(i,j+1) 

A 

S1 (i, j ) 

2  ,  .  ,  % 

2 

S  (1,3+1) 

S  (l,  j) 

+  bu(i, j) 


(IV. 11) 


y(i»j)  = 


R(i , j ) 
S  (i  ,  j  ) 


where : 


C  [b1Q  ...  bnQ  bQ0  0  ...  0  1  0  ...  0] 


b  =  [1  0  ...  0  aQl  ...  aQm  bQ1  ...  bQm]  (input  vector) 


Transition  a  = 
Matrix 


a\m 

:  o 

m 

'  1  i 

o: 

*». 

-*0. 

•  0 

1 

I 

bii 

o  ! 

t 

1 

1 

o 

J>\ m 

b  run 

~bom 

1 

1 

a .  .  = 
ID 

a .  . 
ID 

-  ai0a0j 

b.  .  =  b.  . 
ID  ID 

ai0b0j 

1  <  i 

<  n 

1  <  j  <  m 

1  <  i  <  n 

0  <  j  < 

m 

O 


.  l 

o 


The  expanded  form  of  Eq .  IV- 11  can  now  be  shown  as: 


-V 


»oV 


D.  PROGRAM  AND  EXAMPLES  FOR  ROESSER'S  EQUATIONS  USING 

RUNG'S  MODEL 

This  program  (Appendix  Q) takes  as  initial  conditions  one 
horizontal  state  and  two  vertical  states.  The  order  of  hori¬ 
zontal  states  is  given  by  N  and  the  order  of  the  vertical 
states  by  M. 

We  give  two  examples,  one  for  N  =  2  and  M  =  2  (Example  7) 
(two  orders  for  horizontal  states  and  2  orders  for  vertical 
states)  and  one  for  N  =  4  and  M  =  3  (Example  8)  (four  orders 
for  horizontal  states  and  three  orders  for  vertical  states) . 
The  first  example  is  for  a  matrix  2x2  and  the’ second  example 


1 1 
n 

i 

it 

r> 

e.C'.v'i 

S',  0.0) 

S*-  O*^ 

s\  OS) 

StOA 

s'  0  a') 

s**0^ 

s\G  A 

s\ 

«.iC.vV 

^>0-^ 

<U(?^ 

s\c?,^ 

3tc*A 

S*  Cv*-V 

s;c^ 

s\.(^ 

^  O-r^ 

S\(^ 
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S»C>A^  StC'^X  Si6^\ 
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0-A  «iCs,>\8*(*.^M 

iSt(ki\  SiOA 

S  i  ^  St  SiCjtC\ 

1  1 

^  1  C.v£> ,  £>AA  .*2  fct&Sy  AGaYB^CvS) 
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s»\;.v\  s\^\  sW 

^  i  C\&  „  SAvA  S\CjA 
s\o.^.stC\t)JsWN) 

l 
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S.C\t>,  SxOASiCi^j 

1  s\  Cvb,s\C^\siod 
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S.Ov\  Sv(k^  ,S»(.\A 

SiCi'.A,  StAA  ,SiW 

1 

2-  AA ,  e^A.^c.M.A  .bam^I  bc^a,  BGA 

(s  A  ,  S-.AA. 

i  S*-,  (>sA ,  S\AA , s\(.s^>  O^i),  s\Csi> ,  s\  Cu,^ 

1 

B.  Civ\BtCti  A> ,  Bt^u.^ 

Sji^^i  StC“.^\  SiO*A 
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E.  NUMERICAL  EXAMPLES  FOR  KUNG'S  MODEL 

The  following  presents  three  examples.  The  first  one 
corresponds  to  an  "all-pole"  2-D  low-pass  filter.  The  second 
one  is  an  "all-zero"  2-D  band-pass  filter  (sin  sin  w») • 

The  third  one  is  also  a  band-pass  filter.  All  these  examples 
are  second  order.  The  outputs  of  these  examples  are  produced 
using  Kung's  [Ref.  8]  state-space  model.  In  this  formulation, 
for  a  second  order  system,  we  require  two  horizontal  states — 
Rl(i,j)  and  R2(i,j)  and  four  vertical  states,  Sl(l)(i,j), 

SI (2 )  (i , j ) ,  S2 (1)  (i , j )  and  S2  (2)  (i , j ) .  The  program  listing 
for  implementing  this  model  is  given  in  Appendix  D. 

Example  #9 

The  system  parameters  and  the  initial  conditions  chosen 
for  this  example  are  as  listed  in  Table  4.1.  The  2-D  D.F.T. 
|Y(m,n)|  of  the  output  sequence  y ( i , j )  produced  by  the  program 
in  Appendix  D  is  shown  in  Fig.  4-9a.  The  corresponding  contour 
map  is  shown  in  Fig.  4-9b. 

Example  #10 

The  parameter  coefficients  and  the  initial  conditions  for 
this  example  are  listed  in  Table  4.2.  The  2-D  D.F.T. 
sequence  |Y(m,n) |  for  this  example  is  illustrated  in  Fig.  4-10a 
and  Fig-M.lOb  shows  the  associated  contour  map. 

Example  #11 

The  parameter  coefficients  and  the  initial  conditions  for 
this  example  are  listed  in  Table  4.3.  The  2-D  D.F.T.  sequence 
|Y(m,n)  j  for  this  example  are  illustrated  in  Fig.  4-lla  and 


Figure  4-llb  shows  the  associated  contour  map. 


TABLE  4.1 


NUMBER 

OF  HORIZONTAL 

STATES (N 

®it04) ;  2 

NUMBER 

OF  VERTICAL  STATES (M=ltC4) :  2 

DIMENSION  OF 

OUTPUT ( 1 

toi.3)  :  l 

5 

ENTER 

INITIAL 

CONDITIONS  FOR 

HORIZONTAL  R ( # 

*) 

R  1(1, 

1 )  : 

0 

R  2(1, 

1)  : 

o 

R  1(1, 

2)  : 

0 

R  2(1, 

2)  : 

0 

R  1(1, 

3)  : 

0 

R  2(1, 

3)  : 

0 

R  1(1, 

4)  : 

0 

R  2(  1, 

4)  : 

0 

R  1(1, 

3)  ! 

0 

R  2<  1, 

3)  : 

0 

R  1(1, 

6)  : 

0 

R  2(1, 

6)  : 

0 

R  1(1, 

7)  ; 

0 

R  2(  1, 

7)  : 

0 

* 

R  1(1, 

3)  : 

0 

R  2(1, 

3)  : 

0 

R  1(1, 

9)  : 

0 

R  2d, 

3)  : 

0 

R  1(1, 

10)  : 

0 

R  2(1, 

10)  : 

0 

R  1(1, 

IDs 

0 

R  2(1, 

11): 

0 

R  1(1, 

12)  : 

0 

R  2(  1, 

12)  : 

1.) 

R  Kl, 

13)  : 

0 

R  2(1, 

13)  : 

0 

R  1(1, 

14)  s 

0 

R  2(1, 

14)  : 

0 

R  1(1, 

IS)  : 

0 

R  2(1, 

15)  : 

0 

ENTER 

INITIAL 

CONDITIONS  FOR 

VERTICAL  SI (#. 

#> 

SI  (  1) 

<  1, 

1) 

0 

SI  <  2) 

<  1, 

1) 

0 

SI  (  1) 

<  a. 

1 ) 

0 

SI  <  2) 

<  a. 

1) 

0 

SI  (  1) 

(  a. 

1) 

0 

SI  (  2) 

(  a. 

1 ) 

0 

SI  <  1 ) 

<  4, 

1 ) 

0 

SI  (  2) 

(  4, 

1 ) 

0 

SI  (  1 ) 

<  5, 

1 ) 

0 

SI  (  21 

(  5, 

1) 

0 

SI  (  1) 

(  s. 

1 ) 

0 

SI  <  2) 

<  S, 

1 ) 

0 

SI  (  1 ) 

<  7, 

1 ) 

0 

SI  <  2) 

(  7, 

1 ) 

0 

SI  (  1) 

<  a. 

1 ) 

0 

SI  (  2) 

<  a. 

1 ) 

0 

SI  (  1 ) 

(  3, 

1) 

0 

St  <  2) 

<  3, 

1) 

0 

SI  (  1 ) 

(  10, 

1 ) 

0 

SI  (  2) 

(10, 

1) 

0 

Si  (  1 ) 

(11, 

1 ) 

0 

SI  (  2) 

<11, 

1 ) 

0 

St  '  1  ) 

; ) 

o 
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•j-v.-v-:- 

31 


si  (  o 

SIC  3)  <13,  1)  :  O 
Sic  1)  <  1**,  1)  :  o 
SIC  3)  (  14,  1)  !  o 
SIC  1 )  < IS, l)  :  O 
SIC  3)  (15,  1)  :  0 

ENTER  INITIAL  CONDITIONS  FOR  VERTICAL  S 

sac  i>  c  .i,  i)  s  o 

sac  3)  C  1,1):  0 
sac  DC  3,1):  O 

sac  3)  c  a, i) :  o 

sac  DC  3,  1  >  :  O 
sac  3)  <  3,  1)  :  0 
sac  1)C  4,1):  o 
sac  3)  <  4,  1 )  :  O 
sac  1)C  3, 1) :  0 
sac  3)  C  5,  1)  :  0 
sac  1)C  S,  1 )  :  0 
sac  3)  (  S,  1)  :  0 
sac  DC  7,1):  o 
SEC  3)  C  7,  1)  :  0 

sac  dc  a,  i)  :  o 
sac  a) c  a, l) :  o 
sac  1)C  3,  1)  :  0 
sac  3) C  9, 1) :  o 
sac  1)  ( 10,  1)  :  0 
sac  3)  <10, 1) :  0 
sac  1 )  C 1 1 ,  1)  :  0 
sac  3) <11, D :  O 
sac  D  < la,  d  :  o 

sac  3) < 13,  D :  0 
sac  1)  <13,  D  :  O 
sac  3) <13,  D :  0 
sac  D  <  14,  D  :  0 
sac  3)  <14,  D  :  0 
sac  1>  <15,  1) :  O 
sac  3) (IS, D i  0 

ENTER  VALUES  FOR  THE  INPUT  VECTOR  (#.<*) 

4  <  0  D  :  -O.  33 

a<0  3) i  0 

CKO  1)  :  O 

0(0  3) :  0 


ENTER  ELEMENTS  OF  THE  TRANSITION  MATRIX (#. #> 

a (  10) :  -O. 133 

a<  30):  -0.35 

a  (  1  1)  :  -O.  l 

a  <  3  1)  :  O 

ac  1  3) :  0 

a  <  3  3)  :  -O.  1 

b(  1  1)  :  0 

b<  3  D  :  0 

b<  1  3)  :  0 

b<  3  3)  :  0 

ENTER  VALUES  FOR  THE  OUTPUT  VECTORC#. #> 
b  COO)  :  1 
b (  10)  :  O 
b (  30 ) :  0 
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Figure  4-9a.  2-D  D.F.T.  Sequences,  |Y(m,n) |  for  Example 


Figure  4-9b.  Contour  Map 


TABLE  4.2 


NUMBER  OF  HORIZONTAL  STATES <N= It oh ) : 
NUMBER  OF  VERTICAL  STATES  ( M»  1 1 o1* )  :  2 


DIMENSION  OF  OUTPUT  ( lto25)  :  17 


ENTER  INITIAL  CONDITIONS  FOR  HORIZONTAL  R(t*.  *) 


Rl(l,  1) 

o 

R  3(1,  1) 

0 

R  1(1,  3) 

0 

R  2(1,  2) 

Rl(l,  3) 

0 

R  2(1,  2) 

0 

R  1(1,  4) 

0 

R  2(1,  4) 

0 

R  Id,  5) 

0 

R  2d,  S) 

0 

R  Id,  S) 

0 

R  2  (1,  6) 

0 

r*( 

ft 

o 

R  2  (  1,  7) 

0 

R  1(1,  3) 

0 

R  2 (  1,  3) 

0 

R  Id,  3) 

0 

R  2  d,  3) 

0 

R  Id,  10) 

0 

R  2(1,10) 

0 

R  1  (1,  11) 

0 

R  2(1,11) 

o 

R  1  (1,  12) 

0 

R  2  (1,  12) 

0 

R  1(1,13) 

o 

R  2(1,13) 

0 

R  1(1,14) 

0 

R  2(1,14) 

0 

R  1  (1,  15) 

0 

R  2(1,15) 

0 

Rid,  IS) 

0 

R  2(1,16) 

0 

R  1(1,17) 

0 

R  2(1, 17) 

0 

ENTER  INITIAL  CONDITIONS  FOR  VERTICAL  Sl(#.#> 


SI  ( 

1 ) 

(  1,1): 

0 

SI  ( 

2) 

(  1,1): 

0 

SI  ( 

(. ) 

(  2,  1 )  : 

0 

SI  ( 

_) 

(  2,  1 )  : 

0 

31  ( 

1 ) 

(  3,  1)  : 

o 

SI  ( 

2) 

(  3,  1)  : 

0 

SI  < 

1 ) 

(  4,  1)  : 

0 

SI  < 

2) 

(  4,  1  )  : 

0 

SI  ( 

1 ) 

(  5,  1 )  : 

0 

SI  ( 

2) 

<  5,  1)  : 

0 

SI  < 

1 ) 

(  6,  1 )  : 

<3 

SI  ( 

2) 

(  6,  1 )  : 

0 

SI  ( 

1) 

(  7,  1)  : 

0 

SI  ( 

2) 

(  7,  1 )  : 

O 

SI  < 

l ) 

<  3,  l)  : 

0 

SI  ( 

2) 

(  3,  1 )  : 

0 

SI  < 

1 ) 

(  3,  1)  : 

0 

SI  ( 

2) 

(  3,  1 )  : 

0 

si  ; 

1 ) 

(  10,  1 )  : 

0 

3 1  C 

,£) 

(  10,  1  )  : 

0 

si  ; 

1 ) 

d 1.  1 )  : 

0 
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Figure  4-10b.  Contour  Map  for  Figure  4^10a 


TABLE  4.3 


•MUs'lStt?  OF  1  7  IS  ■'  '-t*  1  5  •  :  1 

VI  X  CAL  3  TP  T. IS  >'  •*»  :•.  :  i. 

OF  OUTPUT  C  ItOiSSJ  :  17 

riPL  CONDITIONS  FOR  HORIZONTAL  RC*. #> 


NUMBER 

OR 

VI** 

D I MENS ION 

QF  ( 

ENTER 

INITIAL 

R  i(l. 

1  )  : 

0 

R  i:d. 

1 )  : 

o 

R  1(1, 

3)  : 

o 

R  3d, 

3)  : 

0 

R  Id, 

3)  : 

0 

R  3  d. 

3)  : 

0 

R  1(1, 

4)  : 

0 

R  3d, 

4)  : 

0 

R  X  d , 

3)  : 

0 

R  3(1, 

3)  : 

0 

R  Id, 

S)  : 

0 

R  3d, 

6)  : 

o 

Rid, 

7)  : 

0 

R  3d, 

7)  : 

0 

R  Id, 

S)  : 

0 

R  3d, 

3)  s 

0 

R  Id, 

3)  : 

0 

R  3d, 

3)  s 

0 

R  Id, 

10)  s 

o 

R  3d, 

10)  : 

o 

Rid, 

U)  : 

0 

R  3d, 

11): 

0 

R  Id, 

.13)  : 

o 

R  3d, 

13)  : 

o 

R  Id, 

13)  : 

o 

R  3d, 

13)  : 

0 

R  Id, 

14)  : 

o 

R  3d, 

14)  i 

0 

R  Id, 

13)  : 

0 

R  3d, 

IS)  ! 

o 

R  1(1, 

16)  : 

0 

R  3d, 

16)  : 

0 

R  1(1, 

17)  : 

0 

R  3d, 

17)  s 

0 

ENTER 

INIT 

I  PL 

SI  (  1 ) 

<  1, 

1 )  : 

SI  <  3) 

<  1, 

I)  s 

SI  <  1 ) 

(  3, 

1)  : 

SI  (  3) 

(  3, 
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Figure  4-lla.  2-D  D.F.T.  Sequences, 


fj  DFT 


Figure  4-llb.  Contour  Map  for  Figure  4-lla 


Once  again,  to  verify  the  correctness  of  our  program. 


the  D.F.T.  |Y(m,n)j  was  compared  to  |h(cj^,u2)|. 

and  the  corresponding  contour  maps  are  shown  in  Fig.  4-12a,b, 

Fig.  4-13a,b  and  Fig.  4-14a,b  for  examples  9,  10  and  11, 

respectively. 

F.  SUMMARY  OF  PROGRAMS  DEVELOPED 

The  programs  which  have  been  written,  cover  the  following 
orders  based  upon  the  different  models. 


Appendix 

Order 

Model 

#_ 

of  States 

A 

1st 

Roesser 

1 

horizontal,  1 

vertical 

C 

2nd 

Roesser 

2 

horizontal,  1 

vertical 

D 

Multi-order 

Rung 

A 

horizontal,  2n 

vertical 

In  order  to  check  the  program  listing,  the  same  first 
order  example  was  used  on  all  programs.  Identical  results 
were  obtained.  Similarly,  identical  second  order  examples 
were  used  in  Programs  C  and  D  and  produced  identical  outputs. 
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Figure  4-13b.  Contour  Map  for  Figure  4-13a 
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Figure  4-14a.  Transfer  Function  jlHz-^z^l, 

for  Example  #11 


Figure  4-14b,  Contour  Map  for  Figure  4-14a 


V.  USE  OP  SSPACK  PACKAGE 


A.  SSPACK 

"SSPACK"  is  a  "state  space  system  package,"  [Ref.  21]  that 
is  an  interactive,  state-of-the-art,  software  package  for  the 
analysis,  design,  and  display  of  one-dimensional  state-space 
systems.  The  work  which  follows  adapts  this  program  so  that 
it  can  be  used  to  produce  2-D  data  fields  from  state  space 
formulations.  A  brief  description  of  SSPACK  fpllows. 

SSPACK  is  useful  for  a  variety  of  applications  in  signal 
processing  and  control  [Ref.  22] .  The  package  consists  of  a 
supervisor  which  controls  the  operation  of  the  software  and  a 
set  of  independent  programs  which  communicate  using  disk  files. 
The  core  of  the  package  are  the  pre-  and  post-processors.  The 
state-space  pre-processor  (SSPREP)  program  aids  in  preparing 
files  for  the  individual  algorithm  programs.  [Refs.  23,24]  Th 
state  space  post-processor  (SSPOST)  program  displays  and  analyz 
the  output  from  the  algorithms.  SSPREP  prompts  with  a  series 
of  questions  in  a  menu  format. 

SSPOST  is  an  interactive  command-drive  processor.  It  is 
designed  to  help  interpret  the  output  of  the  various  SSPACK 
algorithms,  and  display  time  histories: 

A  is  the  Nx  by  Nx  state  transition  matrix; 

B  is  the  Nx  by  Nu  input  transition  matrix; 

C  is  the  Nz  by  Nx  measurement  matrix; 


D  is  the  Nz  by  Nu  feedthru  matrix; 

W  is  the  Nx  by  Nw  process  noise  matrix; 

V  is  the  Nz  by  Nv  measurement  noise  matrix. 

The  SSPACK  works  in  multi-order  form,  using  the  transfer 
function  of  the  1-D  digital  filter. 
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Figure  5-1 


The  present  objective  is  to  use  SSPACK  with  a  2-D  input  data 
field  and  through  the  same  transfer  function,  1-D  digital 
filter,  to  accomplish  2-D  output  data  field. 

B.  DESIGN  OF  2-D  DIGITAL  FILTERS  USING  1-D  DIGITAL  FILTER 

STRUCTURES 

The  idea  of  using  two  types  of  dynamic  elements  is  not  very 
abstract;  it  is  very  natural  in  delay-differential  systems. 
However,  before  considering  its  practical  applications  to 
image  systems,  two  remarks  have  to  be  made.  The  first  is  be¬ 
cause  the  "spatial"  dynamic  elements  seem  unimplementable ,  and 
we  need  to  replace  them  by  time-delay  elements.  Secondly,  in 
order  to  have  a  finite  order,  we  shall  only  consider  a  bounded 
frame  system,  i.e.,  we  assume  that  the  picture  frame  of  interest 


is  an  M  *  N  frame  (with  vertical  width  M  and  horizontal  length 
N) .  Note  that  in  order  to  use  time  delay  elements,  we  need 


first  to  find  a  way  to . code  a  2-D  spatial  system  into  a  1-D 
(discrete  time)  system  and  vice  versa. 

Thus  we  propose  the  following  system,  composed  of  three 
subsystems  in  series: 

i)  The  Input  Scan  Generator  codes  the  2-D  spatial  input 
into  1-D  time  data  according  to  the  mapping  function 

t(-,*)  t ( i , j )  =  iM  +  jN  ,  0  <  i  £  N-l  (V.l) 

0  _<  j  _<  M-l 

where  M  and  N  are  relatively  prime  integers.  for  example,  we 
consider  a  2-D  input  data  u(i,j): 

(0,0)  (0,1) 

(1,0)  (1,1) 

u(i,j)  = 

•  • 

•  • 

•  • 

(N-l , 0 )  (N-l , 1 ) 

Scanning 

The  data  field  u(i,j)  is  scanned  to  produce  u(+)  as  follows 

u  ( i ,  j )  =  (0,0) , (0,1) , (0,2) ,. . . , (0,M-1) , (1,0) , (1,1) , . . . , 

( 1 , M- 1 ) , (N-l, 0) ... (N-l, M-l) 

{ u ( t ) } ,  t  =  0,  1,  2,  M,  M+l,  M+2 ,  ...,  (M-l) (N-l) ,  t  =  iM  +  jN 


For  example, 


I 


u(i,j)  = 


yields 


y  ( t)  =  [1  0  0  0  0  ...  0  0  0  0] 


ii)  A  1-D  (discrete  time)  digital  filter  processes  the 


1-D  data  generated.  This  subsystem  is  implemented  by  replacing 


z. 1  by  6,  z-1  by  A  in  a  2-D  circuit  realization  (e.g.,  2-D 


controller  form).  6  and  A  are  chosen  as: 


6  =  D  =  M-units  delay  element 


A  =  D  =  N-units  delay  element 


iii)  The  Output  Frame  Generator  decodes  the  1-D  (discrete 


time)  output  of  the  1-D  digital  filter  described  above  into 


a  2-D  (discrete-spatial)  picture  according  to  the  inverse 


mapping  of  (V.l) 


(i (t) , j (t) )  =  Pt  Mod  N , [ t- ( Pt  Mod  N)M]/N) 


( V.  2 ) 


where  P  is  a  unique  integer  such  that  PM-PN  =  1  and  0  <  P  <  N. 


This  formula  is  given  in  [Ref.  2],  Alternately,  we  can  compute 


(i / j )  as 


i  =  t  Mod  N 


j  =  Quotient  (t/N) 


» 


?  ji  rw  *.#  ■  •  *  •  ^  .* '  »  • . 


For  example  we  suppose  t  =  19  with  N  =  10  and  M  =  9. 

19 

The  corresponding  value  in  the  2-D  case  will  be  i  =  Remainder 

19 

=  9  and  j  =  Quotient  {jq-}  =  1.  So  in  the  2-D  case  we  will 
have  (i, j)  =  (9,1) . 

Another  Example:  For  M  =  4  and  N  =  5,  the  single  index  •£ 
will  be  mapped  into  (i,j)  as: 


j 


0 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

The  procedure  for  implementing  2-D  filters  using  1-D  filter 
structures  is  as  shown  below  in  Fig.  5-2. 


F\<cLb 


Figure  5-2 


The  index  scanning  is  required  for  the  input  data  so  SSPACK 
can  be  carried  out  simply  because  the  input  is  assumed  to  be 
a  2-D  unit  pulse.  This  is  followed  by  implementing  the 
corresponding  1-D  filter  of  Fig.  5-3  using  SSPACK  to  convert 
the  1-D  output  from  SSPACK  to  a  program  for  output  index 
mapping — written  as  shown  in  Appendix  E.  The  2-D  Fourier 
transform  of  the  resulting  2-D  field  is  then  computed. 

Considering  a  bounded  frame  (MxN)  system  it  is  interesting 
to  know  the  dimension  of  the  global  state  (or  initial  condi¬ 
tions)  needed  to  process  the  MxN  future  data  field.  Since 
vertical  states  convey  information  vertically,  all  the  verti¬ 
cal  states  along  the  X-axis  are  necessary  initial  conditions 
and  their  dimension  is  mN.  Similarly,  all  the  horizontal  states 
along  the  Y-axis  are  necessary  initial  conditions  (with 
dimension  nM) .  They  convey  information  horizontally. 

Therefore,  in  the  bounded  frame  case  a  total  number  of 
mN+nM  are  needed  to  summarize  the  "past"  information.  This 
very  same  idea  can  be  used  again  from  a  computational  point 
of  view.  Indeed,  the  number  of  required  storage  elements  for 
recursive  computations  is  also  equal  to  mN+mN  if  initial  condi¬ 
tions  are  not  zero.  However,  it  is  quite  often  the  case  that 
the  system  starts  with  zero  initial  conditions;  the  size  of 
storage  required  is  reduced  to  mN  (respectively,  nM)  which  is 
used  to  store  the  updated  data  row  by  row  (respectively,  column 
by  column) .  No  storage  is  needed  for  the  rest  of  the  initial 
conditions — nM  horizontal  states  (respectively,  mN  vertical 
states)  since  they  are  assumed  to  be  zero.  This  is  consistent 


with  the  results  of  Read  [Ref.  24]  derived  from  a  direct 
polynomial  approach. 


Another  interesting  observation  concerns  the  dimension 
of  the  1-D  figital  filter  contained  by  our  2-D  digital  filter 
design  discussed  above.  Since  it  needs  nM  unit-delays  and 
mN-unit  delays,  the  corresponding  1-D  state-space  also  has  a 
dimension  equal  to  nM+mN.  Note  that,  despite  the  high  dimen¬ 
sion  of  the  corresponding  1-D  filter,  its  high  sparsity  is  very 
encouraging  for  further  studies.  In  short,  following  the 
above  method  of  designing  a  2-D  filter,  for  the  first  order 
case. 


H  ( zi ' z2 )  = 


1+a10Zl1+a01Z21+allZl  z2 


(V.3) 


Using  the  above  approach  we  get  the  1-D  filter  realization  for 
this  2-D  filter  which  turns  out  to  be  as  shown  in  Fig.  5-3. 

The  detailed  matrix  equations  for  realizing  Eq.  (V.3) 
using  SSPACK  can  be  written  as.  The  SSPACK  produces  a  1-D 
sequence,  which  converted  into  a  2-D  sequence  using  the  output 
index  mapping  formulae  discussed  earlier.  The  listing  of  a 
program  which  does  this  mapping  is  shown  in  Appendix  E. 

After  obtaining  the  valid  2-D  output  data  sequence  y(i,j) 
we  next  compute  its  2-D  D.F.T.  to  produce  |Y(m,n)  j  which  for  this 
example  is  plotted  in  Fig.  5-4a.  The  corresponding  contour 
map  is  as  shown  in  Fig.  5-4b. 
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Fiqure  5-4a.  2-D  D.F.T.  Sequences,  |Y(m,n)|  for  Example 


Figure  5- 4b.  Contour  Mat 


VI.  CONCLUSIONS 


This  thesis  has  dealt  with  the  problem  of  modelling  2-D 
data  fields  in  the  state-space  domain.  First  of  all  we  have 
pointed  out  the  main  problems  associated  with  the  extension  of 

1- D  time-discrete  state-space  models  to  2-D  data  fields.  The 
remaining  part  of  the  thesis  has  been  divided  primarily  in  3 
parts . 

In  the  first  part  we  describe  Roesser's  [Ref.  5]  approach  to 
modelling  2-D  systems  in  the  state  space  domain.  Extensive  com¬ 
puter  simulation  results  are  presented  to  verify  the  functioning 
of  this  approach.  This  modelling  approach  has  been  tried  out 
for  the  scalar  (1  x l)  as  well  as  for  higher  order  (2  x 2)  etc. , 

2- D  systems. 

The  second  part  deals  with  a  modification  of  Roesser's  approach 
as  described  by  Kung  [Ref.  7].  The  main  advantage  of  this  approach 
is  that  the  2-D  state-space  model  can  be  realized  as  a  2-D 
circuit.  More  importantly,  this  2-D  circuit  realization  can  be 
implemented  as  a  1-D  digital  filter.  Computer  simultation  studies 
that  have  been  carried  out  substantiate  the  making  of  this  model. 
The  1-D  filter  realization  obtained  in  this  part  turns  out  to 
be  a  very  convenient  starting  point  for  the  nezt  part  of  our 
effort,  dealing  with  the  use  of  the  1-D  SSPACK  commercial  soft¬ 
ware  package  designed  for  dynamic  system  simulation. 

In  the  final  part  of  the  thesis,  we  make  use  of  the  1-D 
filter  realization  of  2-D  state-space  model  obtained  in  the  second 


part,  and  implement  this  filter  using  SSPACK.  Some  additional 
programming  effort  reuiqred  for  input  and  output  mapping  was 
necessary.  Programs  for  converting  2-D  input  and  output  se¬ 
quences  to  1-D  have  been  written  separately.  In  this 
fashion  we  have  succeeded  in  extending  the  applicability  of 
the  SSPACK  to  simulating  2-D  linear  systems  as  well.  Once 
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THE  PURPOSE  OF  THIS  PROGRAM  IS  TO  COMPUTE  PND  GRAPH  THE 
EQUATIONS  OF  ROBERT  P. ROESSER  IN  THE  "DISCRETE  STATE-SPACE 
MODEL  FOR  LINEAR  IMAGE  PROCESSING". 


EVANGELOS  THEOFILOU 


PROGRAM  2D-DATA-F I ELD 


*****  VARIABLE  DECLARATIONS  ***** 

REAL  R(23.  23)  ,  S  (£5,  23)  ,  R1  (2)  ,  R2(2)  ,  3  (31, 31)  , 

*  RLPART, IMGPART, ZFI31, 31) , VERTEX ( IS)  , ZL2V (31 > 

INTEGER  MASK  ( 3000 )  ,  LD I G  < 3 1 )  ,  LWGT  ( 3 1 ) 

CHARACTER*!  ANSWER 
CHARACTER*20  CTEXT 


DATA 


XLQL/O. 0/, VLOL/0. 0/, XUPR/8.  3/,  YUPR/7.  0/, 
ZLOU/ 1. 0E33/, IPROJ/O/, NRNG/ 100/ 


24  C 


MAIN 


PROGRAM 


**♦♦****•*♦*♦* 


26  C 

*****  ASK  THE 

REQUIRED 

27 

10  WRITE  (*. *)  ’2 

NTER  VAL 

23 

WRITE  (*,333) 

'  Al  :  ’ 

23 

READ  (*, *>  Ai 

30 

WRITS  (*,333) 

•  AS:  ’ 

31 

READ  !*,*>  AS 

32 

WRITE  (*,333) 

'A3:  ’ 

33 

READ  (*, *1  A3 

34 

WRITE  (*,333) 

*  A4;  * 

33 

READ  (*, *)  A4 

36 

WRITE  (*,333) 

’91:  ’ 

37 

READ  <  *,  *)  B 1 

33 

WRITE  (*,333) 

’  B£ :  ’ 

33 

READ  (*, *)  B2 

40 

WRITE  (*,333) 

’  Cl  :  ' 

41 

READ  (*, *)  Cl 

42 

WRITE  (*,333) 

’  C2:  ’ 

42 

READ  (*, *)  02 

44 

3  WRITE  (*,402) 

43 

READ  (*, *)  N 

46 

IF  (N  .  GT.  25 

GOTO  3 

47 

48 

WRITE  (*,211) 

’ ENTER  ’ 

DO  33 
WRITE 
READ 


I  *  1,  N 

(*,  403)  ’R(l, 

(*,  *)  R  <  1 ,  I) 


93  CONTINUE 


WRITE  (*, 21 1 )  ’ENTER  *  ,  N,  ’ 
DO  100  I  »  1 , N 
WRITE  (*,404)  '  S(,  ’  ,  I,  ’  1)  : 

READ  (*. *)  S( I, 1 ) 

100  CONTINUE 


INITIAL  CONDITIONS  FOR  MATRIX  S(».#)' 
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65  C 
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2  70 

2  71 

2  72 

2  73 

2  74 

2  75 

76 

77  C 

78 

1  73 

2  80 

2  81 

82 

83  C 

84 

1  85 

2  86 

2  87 

88 

33  C 

30 

31 

32 

1  33 

1  34 

1  35 

36 

37 

38 

33 
100 
101 

102  C 

103 

104 

105 

106 

107 

108 
103 
1 10 
111 
112 

113 

114 

115 
1  IS 
117 
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WRITE  (*,413) 

READ  (*,200)  ANSWER 

IF  ((ANSWER  .EQ.  ’V’)  .OR.  (ANSWER  .ED.  'y’>>  GOTO  10 

U  =  1.0 

*****  COMPUTE  R  PND  S  MATRICES  ***** 

DO  101  I  =  1,N 

DO  101  J  =  1,N 

IF  (1  +  1  .LE.  N)  THEN 

R ( 1*1, J>  =  P1*R(I,J)  *  P£*S ( I , J>  +  8 1 *U 
END  IF 

IF  (J+l  .LE.  N)  THEN 

S(I,J*1)  =  P3*R ( I, J)  +  A4*S(I,J>  +  B£*U 
END  IF 
U  =  0.  0 

101  CONTINUE 

****  FILL  0’ s  THE  TWO  DIMENT IONPL  GRID  OF  CONTROL  POINTS  **** 
DO  102  I  =  1,31 
DO  102  J  =  1,31 
Z(I,J)  =0.0 

102  CONTINUE 

*****  COMPUTE  Z  MATRIX  ***** 

DO  103  I  =  1 , N 
DO  103  J  =  1 , N 

Z(I,J>  =  C1*R(J,  I)  *  C2*S  <  J,  I ) 

103  CONTINUE 

*****  OUTPUT  THE  Z  MATRIX  ***** 

WRITE  (*,205)  ******  Z  MATRIX  ’  ,  N, ’  X  *  ,  N, '  ****** 

WRITE  (*,212) 

DO  104  I  =  1 , N 

WRITE  (*,300)  (Z(I,J>,  J  =  1 , N) 

WRITE  (*,210) 

104  CONTINUE 
WRITE  (*,213) 

WRITE  (*, 418) 

READ  (*,200)  PNSWER 

IF  ((ANSWER  . NE.  ’ Y* >  .AND.  (PNSWER  . NE.  ’ y’ ) )  GOTO  18 

*****  ASK  THE  PARAMETERS  FOR  THE  GRAPH  ***** 

15  WRITE  (*,210) 

WRITE  (*, *)  ’ ***  ENTER  PLOT  PARAMETERS  *** 
WRITE  (*,405) 

READ  (*, *)  AZIM 

WRITE  (*,406) 

READ  (*, *)  ELEV 

WRITE  (*,408) 

READ  (*, *)  ITRIM 

WRITE  (*,403) 

READ  (*, *)  IDIV 

WRITE  (*,411) 

READ  (*, 133)  CTEXT 
WRITE  (*,401) 

READ  (*,200)  ANSWER 

*****  INITIALIZE  PLOT 33  ***** 
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151  C 
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160 
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IF  ((ANSWER  .  EQ.  ’Y*>  -OR.  (ANSWER  .Ed.  ’  y’  )  >  THEN 
CALL  PLOTS <0,0, 2) 

ELSE 

CALL  PLOTS (0, 33, 33) 

END  IF 

CALL  WINDOW (XLOL, YLOL, XUPR, YUPR) 

*****  DRAW  THE  MESH  SURFACE  OF  THE  GRAPH  ***** 

CALL  ME3HS  ( Z,  31, 31,  N,  N,  AZIM,  ELSV,  0.  3,  0.  3,  7.  5,  3.  3,  ID  IV,  0, 

*  3, IPROJ, 1, ZLOW, 3, ITRIM, MASK, VERTEX) 

*****  ANNOTATION  OF  THE  GRAPH  ***** 

CALL  SYMBOL (3. 5, O. 3, 0. 2, ’ AZIMUTH:  *,0.0,10) 

CALL  NUMBER (333. 0,  333.  0,  0.  2,  AZIM,  O.  0,  2) 

CALL  SYMBOL (5. 5, 0. 0, 0. 2, *  ELEVATION : ’  , 0. 0,  10) 

CALL  NUMBER (333. O, 393.  0, 0.  2,  ELEV, O. 0,  2) 

DY  «  <Z<1, 1 1/90.0)  *  ELSV 

CALL  P3D2D ( 1. 0,  1.0, Z  <1,  1) -DY, XR, YR) 

CALL  SYMBOL(XR, YR, 0.  23,  ’  ,  0.  0,  1) 

CALL  SYMBOL < 1. 0, 0.  1, O. 2,  ’ *  -  ORIGIN’ , 0. 0,  10) 

CALL  SYMBOL ( 1. O, 6. 73,  0.  25, CTEXT, 0.  0,  20) 

CALL  SYMBOL (6. 0, 6. 3, 0. 2, ’ 2-D  DATA  FIELD’ , 0. 0, 14) 

*****  OUTPUT  THE  GRAPH  ***** 

CALL  PLOT  (0.  0,  0.  0,  333) 

WRITE  (*,412) 

READ  (*,200)  ANSWER 

IF  ((ANSWER  .  EQ.  '  Y*  >  .OR.  (ANSWER  . EQ.  ’ y’  ) >  GOTO  IS 

13  WRITE (*, 417) 

READ (*,200)  ANSWER 

IF  ((ANSWER  .EQ.  ’  Y’  >  .OR.  (ANSWER  . EQ.  ’ y’  ) )  THEN 

****  FILL  O’ s  THE  TWO  DIMENTIONAL  GRID  OF  CONTROL  POINTS  **** 

DO  106  I  »  1,31 
DQ  106  J  =  1,31 
ZF ( I , J)  =0.0 
106  CONTINUE 

ZFMAX  -  -3.  3E20 
ZFMIN  =  3. 3E20 
DN  »  (N— 1 ) /2. O 
P  =  6. 2S3185 
DO  107  I  *  1,N 
DO  107  J  =  1 , N 
RLPART  »  O. O 
IMGPART  =  0.0 
DO  108  L  =  1 , N 
DO  108  K  =  1,N 

R 1  ( 1 )  =  COS ( — P* (L— 1 >  * ( I — DN— 1 ) /N) 

R1  (2)  =  SIN (— P*  <L— 1 ) * ( I— DN— 1 > /N) 

R2 ( 1 )  =  COS (-P* (K-l ) *( J-DN-1 ) /N> 

R2 (2)  =  SIN ( — P*  < K— 1 ) *  <  J-DN— 1 ) /N) 

RLPART  =  RLPART  +  Z (L, K) * < R1 ( 1 ) *R£ ( 1 ) 

*  — R 1  (2)  *R2  (2)  ) 

IMGPART  =  IMGPART  +  Z(L, K)*(R1(1) *R2 (2) 

*  +R1 (2) *R2 ( 1) ) 

108  CONTINUE 

ZF ( I , J )  =  SCRT (RLPART**2  *  IMGP0RT**2> 

IF  ( ZF ( I, J)  .  GT.  ZF^ex:  “HEN 
ZFMAX  =  ZF'I.J) 


125 
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^in*»  1 

7 

Microsoft 

FQRTRAN77  V3. 

£ 

178 

END  IF 

a 

173 

IF  ( ZF ( I , J)  .LT.  ZFMIN) 

THEN 

180 

ZFMIN  =  ZF ( I , J> 

181 

END  IF 

* 

182 

107 

CONTINUE 

183 

184  C 

185 

WRITE  (*,205)  '  ***  FOURIER  TRANSFORMATION  > 

y 

X 

2 

186 

WRITE  (*,212) 

187 

DO  103  I  ■  1,N 

i 

188 

WRITE  (*,300)  CZF(I,J),  J 

»  1,  N) 

i 

183 

WRITE  (*,210) 

i 

130 

103 

CONTINUE 

131 

WRITE  (*,213) 

132 

133 

WRITE  (*,418) 

134 

READ  (*,200)  ANSWER 

135 

IF  ((ANSWER  . NE.  1  Y’  )  .AND. 

(ANSWER  . NE.  '  y' 

) >  GOTO  16 

136 

a;e 

Sl-ZZ- 85 
12:34:27 
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137  C 
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THE 
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1 38  30 

WRITE 

<*, 

210) 

133 

WRITE 

(*, 

*)  ’ 

***  E  N 

200 

WRITE 

<*, 

405) 

201 

READ 

<*, 

*) 

AZIM 

202 

WRITE 

<*, 

406) 

203 

READ 

(*, 

*) 

ELEV 

204 

WRITE 

(*, 

*08) 

205 

READ 

<*, 

») 

ITRIM 

206 

WRITE 

(*, 

403) 

207 

READ 

(*, 

* ) 

ID  IV 

208 

WRITE 

(*, 

411) 

203 

READ 

(*, 

133) 

CTEXT 

210 

WRITE 

(*, 

401 ) 

211 

READ 

(*, 

200 ) 

ANSWER 

Si 2  C 


*****  INITIALIZE  PL0T88  ***** 


214 

213 

216 

217 

218 
213 
220 
221 
222 

223 

224 


IF  ( (ANSWER  . ED.  ’ Y’ ) 
CALL  PLOTS  CO, 0,  2) 
ELSE 

CALL  PLOTS  < 0, 33, 33 ) 
END  IF 


,  OR.  (ANSWER  . EQ.  ’ /’  !  )  THEN 


WRITE  (*,420) 

READ  (*,200)  ANSWER 


CALL  WINDOW ( XLOL,  YLOL,  XUPR,  YUPR) 


226 

227 

228 
223 

230 

231 

2*7  0 


1  10 


IF  ((ANSWER  .ED.  ’  Y’  >  .OR.  (ANSWER  . EQ.  ' 
DLEV  =  (  ZFMAX  — ZF-  1 1N)  /FLOAT  ( N ) 

CALL  ZLEVEL (ZF, 31, 31,  N,  N.  DLEV,  ZLEV, N+ 1 > 
DO  110  I  -  1,N+1 
LDIG(I)  -  2 
LWGT(I)  =  l 
CONTINUE 


)  ) 


THEN 


Z33  C 


234 
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236 


*****  draw  THE  CONTOUR  NAP  ***** 

CALL  ZCNTUR  (  ZF,  31 , 3 1 ,  N.  N,  0.  5,  0.  5,  7.  3,  5. 

N*l, 0. 10, 10) 

CALL  SYMBOL  5. 5, 0. o, o. 2. ’ CONTOUR  'IBS’  ,  0.  1 1  ’ 


LEV.  LDIG.  L^GT, 
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ELSE 

*****  DRAW  THE  MESH  SURFACE  OF  THE  GRAPH  ***** 

CALL  MESHS (ZF,  31, 31,  N,  N,  AZIM, ELEV, 0. 5, 0. 5, 7. 5, 5. 5,  IDIV 
3,  I  PRO J,  1, ZLCW, 3,  ITRIM, MASK, VERTEX) 

*****  ANNOTATION  OF  THE  GRAPH  ***** 

CALL  SYMBOL (5. 5,  0.  3,  0.  £,’ AZIMUTH:  >,0.0,10) 

CALL  NUMBER <333. 0, 333. O,  0. £, AZIM, O. 0, £) 

CALL  SYMBOL  <5. 5, 0. 0, O. £,  ’ ELEVATION: > , 0. 0,  10) 

CALL  NUMBER  <333.  O,  333.  0,  O.  £,  ELEV,  0.  0,  £) 

DY  =•  <ZF  ( 1,  1 )  /30.  0)  *  ELEV 

CALL  P3D£D< 1. 0,  1. 0, ZF  < 1, 1) -DY, XR, YR) 

CALL  SYMBOL  <XR, YR, 0. £5, *  *'  , 0. 0,  1 ) 

CALL  SYMBOL <1.0, 0. 1,0. £, >*  =  ORIGIN’ , 0. 0, 10) 


£50 

END  IF 

£51 

CALL  SYMBOL  <1.0, 6. 73, 0.  £5,  CTEXT, 0. 0, £0) 

£53 

CALL  SYMBOL <6.0, 6.5,0. £,  ’  £-D  DFT’,0.0,7) 

£54 

£55  C 

*****  OUTPUT  THE  GRAPH  ***** 

£56 

CALL  PLOT  <0. 0, 0.  0,  333) 

£57 

WRITE  <*, 412) 

£58 

READ  <*, £00)  ANSWER 

£53 

IF  ((ANSWER  .EQ.  ’  Y’  >  .OR.  (ANSWER  . EQ. 

y>  >  ) 

GOTO  30 

£60 

16 

END  IF 

£61 

WRITE  <»,413> 

£6£ 

READ  <  *, 200)  ANSWER 

£63 

IF  ((ANSWER  .EQ.  > Y’  )  .OR.  (ANSWER  . EQ.  ’y 

)  )  GOTO 

10 

£64 

STOP 

£63 

£66 

133 

FORMAT (A£0> 

£67 

200 

FORMAT < A) 

£68 

£05 

FORMAT </, £0X, AES,  I£,  A3,  I£,  A8, /) 

£63 

210 

FORMAT (/) 

£70 

£1 1 

FORMAT </, AS,  12,  A47) 

£71 

£l£ 

FORMAT  </, £X,  ’  (AZIMUTH  3£ 0. 0) ’  , 46X,  ’  (AZIMUTH 

0)’, 

/> 

£73 

213 

FORMAT ( / , EX,  ’  (AZIMUTH  050. 0 ) ’  , 46X ,  >  (AZIMUTH 

140. 

0)’  , 

/) 

£73 

300 

FORMAT ( 10 <F7. 2, IX) ) 

£74 

333 

FORMAT </, 5X,  A4,  \) 

£73 

400 

FORMAT <3X, \> 

£76 

401 

FORMAT </, 3X,  ’ SEND  GRAPH  TO  THE  PR I NTER ( Y  or 

N)  ; 

’  ,  \> 

£77 

402 

FORMAT (/, 3X, ’ NUMBER  OF  ROWS/COLUMNS  FOR  R  AND  S(1  to 

£5) 

£73 

403 

FORMAT  <SX,  A4,  12, A3,  \ ) 

373 

404 

FORMAT  <5X, A2,  12, A5, \ ) 

£80 

403 

FORMAT </, 5X, ’ AZIMUTH <0.  0  to  360.0  DEGREES): 

’ ,  \> 

£81 

406 

FORMAT  < / , 3X , ’ ELEVATION  <  30. 0  to  -30.0  DEGREES):  ’ 

,  ') 

£8£ 

408 

FORMAT  </,  3X,  '  TRIM (O— NO,  l  =  <s,£=Ys):  ’,\) 

£83 

403 

FORMAT < /, 3X,  ’  2,  4  OR  8  SUBGRIDS:  ’,\> 

£84 

41  1 

FORMATt/, 5X,  ’ TITLE  OF  GRAPH(UP  TO  20  CHAR): 

*  ,  \) 

£35 

41£ 

FORMAT </, 3X, > DO  YOU  WANT  TO  CHANGE  PARAME7E 

R3 * 

, 

£86 

413 

FORMAT  </,  3X,  ’  DO  YOU  WANT  TO  REPEAT  THE  PROCESS'1 

£87 

417 

FORMAT </, 5X,  ’ DO  YOU  WANT  FOURIER  TRANSFORMATION 

* 

* 

\  > 

£88 

4ia 

FORMAT </, sx,  ’ DO  YOU  WANT  TO  MAKE  GRAPH  ?  ’, 

\> 

£83 

413 

FORMAT  (/,  3X,  ’  DO  YOU  WANT  TO  CORRECT  t  ’  ,  \  ) 

£30 

4£0 

FORMAT  <  /  ,  5X,  ’  DO  YOU  WANT  CONTOUR  MAP  "•  >,\) 

£31 

END 

Name  Type  Offset  P  Class 

A 1  REAL  £6 


127 


D  Lire#  1  7 

05  REAL 

30 

i  c. :  :  d  7 

Microsoft  P0RTR0N77  V3.  20  02/ Q^* 

03 

REAL 

34 

04 

REAL 

3S 

ANSWER 

CH0R*1 

74 

QZIM 

REAL 

114 

B 1 

REAL 

42 

62 

PEAL 

46 

Cl 

real 

30 

> 

C2 

REAL 

34 

cas 

CTEXT 

CHAR*£0 

126 

intrinsic 

DLEV 

REAL 

216 

DN 

real 

166 

DY 

REAL 

146 

ELEV 

floot 

REAL 

1 13 

intrinsic 

I 

INTEGER*^ 

60 

ID  I V 

INTEGER*^ 

124 

IMGPOR 

REAL 

130 

I  PRO  J 

INTEGER*£ 

22 

ITRIM 

INTEGER*^ 

122 

J 

INTEGERS 

66 

K 

INTEGER*^ 

SOS 

L 

INTEGER*2 

134 

» 

LDIG 

INTEGERS 

6000 

LARGE 

LUlGT 

INTEGERS 

6062 

LARGE 

MASK 

INTEGER*^ 

0 

LARGE 

N 

integer*^ 

58 

NRNG 

INTEGER*^ 

24 

p 

REAL 

170 

R 

real 

o 

LARGE 

91 

REAL 

0 

LORGE 

R2 

REAL 

a 

LARGE 

RLPORT 

REAL 

186 

s 

REAL 

2300 

LARGE 

SIN 

INTRINSIC 

SORT 

U 

real 

76 

INTRINSIC 

vertex 

real 

0 

LORGE 

XLOL 

real 

2 

XR 

REAL 

150 

XUP9 

REAL 

1C 

YUCL 

real 

6 

YR 

REAL 

154 

yupr 

REAL 

14 

Z 

real 

5000 

LORGi- 

ZF 

REAL 

6844 

LORGE 

ZFMAX 

REAL 

158 

ZFMIN 

REAL 

162 

ZLEV 

REAL 

12688 

LPRGE 

ZLOW 

REAL 

18 

Name 

Tyne 

S  i  :s 

Class 

MAIN 

PROGRAM 

MESAS 

SUB9CUT INF 

NUMBER 

SUBROUTINE 

128 

7 


D  Line#  1 

P3D2D 

PLOT 

PLOTS 

SYMBOL 

WINDOW 

ZCNTUR 

Z LEVEL 


Pass  One 


P> 


SUBROUTINE 

SUBROUTINE 

SUBROUTINE 

SUBROUTINE 

SUBROUTINE 

subroutine 

SUBROUTINE 


No  Errors  Detected 
231  Source  Lines 


Microsoft  PQRTRPN77 


Pace  I 
03-EE-as 
1 2 : 34 : 27 
V3.  EG  02/ Q; 


f •}  fu  fu  *“*  f'j  iw  iu 


appendix  b 


0  Lina*  I  7 

1  ^STORAGE:  2 
£  SPAGESi 22:53 


-QRTRAN77  v; 


-aq»  : 

.;:g— i 

21:  i  j 

■  20  55/9* 


3> 

4 

c 

c 

* 

**************************************************** 

* 

3 

r 

*  THE  PURPOSE  OF  "HIS  PROGRAM  IS  TO  COMPUTE  PND  CRASH  "HE 

♦ 

6 

c 

*  FREQUENCY  RESPONSE  OF  A  2-D  DIGITAL  FILTER. 

* 

7 

U 

♦ 

♦ 

3 

C 

* 

evangelos  thecfilou 

* 

3 

c 

10 

c 

PROGRAM  2D— DATA— FIELD 

11 

12 

c 

*****  VARIABLE  DECLARATIONS  ***** 

13 

REAL 

A  (7,  7)  ,  3  (7,  7)  ,  Rl  <7,  7,  2)  ,  R£(7,  7,  2)  , 

14 

* 

rlpart, IMGPART, Z (31, 51) , 

IS 

♦ 

VERTEX ( 16) , ZLEV (31) 

16 

INTEGER 

MASK ( 3000 ) , L3 I G ( 3 1 ) , LWGT (31) 

17 

CHARACTER*  l 

answer 

13 

CHARACTER**' 

0  CTEXT 

13 

* 

20 

DATA 

XLOL/0. 0/, YLCL/O.  0/,  XUPR/3. 5/, YUPR/7. 0/, 

3  1 

* 

ZLCW/ 1.  0E23X ,  IPRQJ/O/, NRNG/ 100/ 

3? 

c 

**  main  program  ************* 

24 

z-m 

10  WRI"E  (*,401) 

26 

READ  <  »,  *> 

IT 

27 

IF  (IT  .3T. 

23)  GOTO  10 

29 

HRITc  (*,  40 

2) 

23 

READ  ( *, ♦ ' 

K 

20 

K  *  K  *  1 

31 

WRITE ( *, *) 

’  ENTER  VALLES  OF  CCEr F I C I ENTS : ’ 

22 

DO  100  I  * 

0,  K-l 

ww 

DO  100  J 

=  0, K-l 

24 

WRITE i * 

■.404)  ’90  ,  I,  ’  ,  ’  ,2,  ’  )  :  ’ 

26 

27 

a 

23 


READ  (*,*)  8dn,W) 
100  CONTINUE 

DO  101  I  »  0,  K-l 
DQ  101  J  =  O,  K-l 


40 

UIRI~c(*,  404)  ’  A  ( •  ,  I,  ’  ,  ’  ,  J,  ’  )  :  ’ 

41 

READ  (*,  *)  A (1*1,  J*1 ) 

4£  101 

CONTINUE 

44 

WRITE  (*,413) 

43 

READ  (*,200)  ANSWER 

46 

IF  ((ANSWER  . EG.  ’ Y’  )  .OR.  (ANSWER  . EQ.  ’y’>>  GOT 

o 

o 

46  C 

****  FILL  O’ 3  THE  TWO  DIMENTIONQL  GRID  OF  CONTROL 

POINTS 

49 

DO  107  I  =*  1 ,  S 1 

30 

DO  107  J  *  1,31 

31 

Z(I,J)  =0.0 

130 


u3  O'* 


34 

53 

36 

57 

38 

38 

50 

51 

52 

53 

54 
65 
So 

67 

68 
S3 

70 

71 

72 

73 

74 

75 

76 

77 
73 
73 
80 
81 
82 

33 

34 
83 
86 

87 

88 
83 

30 

31 

32 

33 

34 
33 

36 

37 

38 
33 

1  oo  c 

101 

102 


107  CONTINUE 

ZMIN  =  3.  3E20 
ZMAX  =  -9. 3E20 
0  =  3. 14133 

STEP  =  2*0  /  (IT-1) 

W1  -  -P  -  STEP 
L  =0 

DO  102  I  =»  1,  IT 
W1  a  ui  +  STEP 
W2  =»  -P  —  STEP 
DO  103  J  a  1, IT 
L  =*  L  +  1 
U2  =  W2  +  STEP 
DO  104  M  a  0,  K-l 


Microsoft  PCRTP0N77  VC. 20  02/34 


1  104  N  a  0, 
R1  (M+l,  N+l, 

K-l 

1 )  a 

COS ( — M 

* 

W1  ) 

R1 (M+l, N+l, 

2)  a 

SIN  (-M 

* 

W1  ) 

R£ (M+l,  N+l, 

1  )  = 

COS (-N 

* 

W£> 

R2 (M+l, N+l, 

2)  a 

SIN ( — N 

* 

W£) 

CONTINUE 
RLNOM  =  O.  O 
IMGNOM  a  0.0 
RLDEN  «•  0. 0 
IMGDEN  =  O. 0 
DO  105  M  =  O 


0.  0 
0.  0 
O.  0 

=  0, K-l 


DO  105  N  =  O,  K-l 
RLNOM  a  RLNQM 


IMGNOM 

RLDEN 

IMGDEN 

CONTINUE 
ELEMENT  a  S 


RLNQM+8  l  M+l,  N+l  )  *(R1  <M+1,  N+l,  1  )  *R2  (M+l  ,  N+l,  1  ) 

-  R1  (M+l,  N+l ,  2) *R2 (M+l , N+l , 2) ) 
I MGNQM+8  (M*  1 ,  N+l )  *  1 R1  (M+l ,  N+l ,  1 )  *R£.(M+1 ,  N+l , 

*  R£  (M+l,  N+l,  1 )  *R1  (M-i-1,  N+l , 
RLDEN+CXM+l,  N+l) *(R1 (M+l,  N+l,  1 ) *R£ (M+l,  N+l,  1 

-  R1 (M+l, N+l, 2) *R£ (M+l,  N+l,  2 
IMGDEN+A (M+l,  N+l) * ( R1  (M+l ,  N+l,  1 ) *R2 (M+l,  N+l, 

+  R2 (M+l, N+l, 1 ) *R1 (M+ 1 , N+l , 


ELEMENT  a  SORT ( RLN0M++2  +  IMGNCM**2)  / 

*  SORT < RLDEN**2  +  IMGDEN**2) 

Z(I,J)  »  ELEMENT 
IF  <  Z ( I , J)  .GT.  ZMOX)  THEN 
ZMOX  a  Z  (  I,  J) 

END  IF 

IF  <  Z ( I , J)  .UT.  ZMIN)  THEN 
ZMIN  a  z ( I, J) 

END  IF 

103  CONTINUE 
102  CONTINUE 

*****  OUTPUT  THE  Z  MATRIX  ***** 

WRITE  (*,203)  ’**+**  Z  M  «  T  R  I  X  ’  ,  IT,  •  X  •  ,  IT,  *  *+***< 

WRITE  (*,212) 


fa  ro  ^  ~  r».i  ft.* 


1 

1 


1,  IT) 


.  r<e** 
103 
10* 

105 

106 
107 

ioa 

109 

1 1 0 
111 
112 

113 

114 


□0  106  I  =  1, IT 

WRITE  (*,300)  (Z(  I,-),  J  = 

WRITE  (*,310) 

106  CONTINUE 

WRITE  (*,313) 

WRITE  (*,416) 

READ  (*,300)  ANSWER 

IF  ((ANSWER  .NE.  '  V’  >  .AND.  (ANSWER  . NE.  ’  y>  )  )  GOTO 


*****  ASK  THE  PARAMETERS  FOR  THE  GRAPH  ***** 


115  20 

WRITE 

210) 

ue 

WRITE 

<*, 

*>  1 

****  E  N 

1 17 

WRITE 

410) 

116 

READ 

<*, 

*) 

AZIM 

1  19 

WRITE 

(  *, 

411) 

120 

READ 

(*, 

*) 

ELEV 

121 

WRITE 

<*, 

413) 

122 

READ 

(*, 

* ) 

ITRIM 

123 

WRITE 

(*, 

414) 

124 

READ 

(*, 

*> 

ID  IV 

125 

WRITE 

<*, 

415) 

126 

READ 

(*, 

193) 

CTEXT 

127 

WRITE 

(*, 

451) 

126 

READ 

(*, 

200  > 

ANSWER 

129 

l-J 

G 

n 

***** 

INITIALIZE  PL0T36 

131 

IF  < (ANSWER 

.  EQ.  ’  Y’  ) 

133 

133 

134 

135 

136 

137 

138 

139 


PARAMETER 


call  PLOTS (0, 0,  3) 
ELSE 

CALL  PLOTS (0,99,  99) 
END  IF 

WRITE  (*,430) 

READ  (*,300)  ANSWER 


l 

1 

1 


140 

141 

142 

143 

144 

145 

146 

147 
146 

149  C 

150 

151 

152 

153 


CALL  WINDOW (XLOL, YLOL, XUPR, YUPR) 

IF  ((ANSWER  . EQ.  ’ Y» >  .OR.  (ANSWER  .ED.  ’ y’ ) >  THEN 
DLEV  =  (ZMAX-ZMIN) /FLOAT ( IT) 

CALL  ZLEVEL  <  Z, 31,51,  IT,  IT, DLEV,  ZLEV,  IT+1 ) 

DO  106  I  =  1,  IT*1 
LDIG ( I )  =2 
LWGT < I )  =1 

106  CONTINUE 

*****  DRAW  THE  CONTOUR  MAP  ***** 

CALL  ZCNTUR ( Z, 51,  51,  IT,  IT,  0.  5,  0.  5,  3.  35,  6.  5,  ZLEV, LDIG, LWGT, 
*  IT+1, 0. 10, 10) 

CALL  SYMBOL (5. 5, O. 0,  O. 2, ’ CONTOUR  MAP’  , 0. 0,  1 1 ) 

ELSE 


vr 

** 

& 

£ 


p* 

.>3 


Is 


■v** 

os 

&• 


PK 


.  r ie-»? 

:sa 
1 


*****  DRAW 


cc 

mJ 

CALL 

MEShS ( 2,  51, 51,  IT,  IT, , 

5£ 

* 

3,  I  PRO J,  1. ZLOW 

-H  /  W 

■  ANNOTATION  OF  THE  G! 

53 

CALL 

SYMBOL  (5.  5,  0.  3,  0,  0,  ’ 

53 

CALL 

NUMBER (333.  0, 333.  0,  0, 

60 

CALL 

SYMBOL  (5.  5,  0.  0,  0.  0,  ’ 1 

£  J 

CALL 

NUMBER (393. 0, 393.  0, 0, 

60 

DY  = 

(  Z  ( 1 ,  1 ) / 30. 0 )  *  ELEV 

63 

CALL 

P3DED  <  1 . 0,  1. 0,  Z  <  1,  1  > 

6* 

CALL 

SYMBOL <  XR,  YR,  0. 05,  ’ * 

65 

CAL_ 

SYMBOL  d.O,  0.  1,  0.  0,  • 

Microsoft  rCn~3h.\77  V3.  ,-0 
HE  m£3h  SURFACE  OF  THE  GRAPH  ***** 

*1,  ELEV,  0.  5,  0.  5,  8.  05,  5.  5,  IDIV,  0 


***** 


.  0,  0) 

:  ’  ,  0. 

.0,0) 


0.  0,  10) 


O,  10) 


END  IF 


ORIGIN’ , 0. 0, 10) 


oa 

S3 


175 


1  7  9 
1  7R 


15 


C3Ui_  SYMBOL  (1.0,  6.  75,  0.  05,  CTSXT,  O.  0,  00) 

COI —  SYMBOL (6. 0, 6. 5, 0. 0, ’ 3-D  DATA  FIELD’ ,0.0,  14) 

*****  OUTPUT  THE  GRAPH  ***** 

C3Lw  PLOT (0. O,  0.  0,  393) 

HR  I  TG  (  *,  *  1 6 ) 

READ  (*,000)  ANSWER 

if  ((Answer  .ec.  • y’ >  .or.  (answer  . eq.  ’y’>>  gc~o  oo 

WRITE  (*,417) 

READ  !*, 000>  ANSWER 

IF  ( (ANSWER  . EQ.  ’ Y’ )  .OR.  (ANSWER  . GO.  ’y’>>  30~0  10 


1  30 
q  < 

S’OP 

*  LJ  * 

130 

199 

-CRmAT  <030  > 

133 

zoo 

FORMAT (A) 

13* 

Z  05 

CQRMAT ( /, OOX, AOS,  10, A3,  10, A8, / ) 

185 

210 

FORMAT < ) 

1 86 

2:  i 

FORMAT (/,  5X,  A60) 

137 

212 

FORMAT  (/,  OX,  '  (AZIMIJ-H  300.  0)  ’  ,  *6X,  ’  (AZIMUT 

H  220.0> 

’  f  /  ) 

iaa 

2 1  7 

FORMAT ( /, OX, ’ ( AZIMO-H  050. 0) ’ , *SX, ’ (AZIMUT 

H  1 40.  0 ) 

’ ,  /) 

183 

300 

FORMAT  < 10 (F7. 0,  IX)) 

130 

*►00 

FORMAT (3X, \ ) 

131 

*♦5 1 

FORMAT  (/,  =x,  ’  SEND  GRAPH  TO  ~HE  PRINTER  (Y 

r  V)  z  • 

\ ) 

130 

401 

F0RMA7</,  5X.  ’ DIMENSION  OF  OUTPUT  MATRIX (1 

to  25)  : 

’ ,  \ ) 

133 

402 

FORMAT*/,  5X,  *  ORDER  OF  TRANSFER  FUNCTION (0 

to  4>  :  1 

1  \J 

13* 

404 

FORMAT  <5X, AO,  I 1, A,  I 1, A3,  \ ) 

1  95 

4  1  0 

FORMAT  (/,  sx,  ’  AZ  IMUTH  <0.  ,;>  to  360.0  DEGREES) 

:  ’  ,  \  > 

136 

41 : 

FCPMAT </,  5X, ’ ELEVATION (30.  0  to  -30.0  DEGRE 

25) :  * , \ 

) 

137 

413 

FORMAT  (/,  5X,  ’  TRIM(  0=*N0,  1  =  <  s,  6=  Vs  )  :  ’  ,  \ ) 

138 

4  14 

FORMAT (/, 5X,  ’ E,  *  OR  3  SUBGRIDS:  ',\> 

133 

415 

FORMAT (/, 5X, ’ TITLE  OF  GRAPH (UP  TQ  GO  CHAR) 

:  ’  ,  \> 

000 

4  1  & 

FORMAT (/, 5X, ’ DO  YOU  WANT  TO  CHANGE  PARAME” 

E.9S  **  ’  , 

\ ) 

00  1 

417 

FORMAT  (/,  SX,  ’  DO  YOU  WANT  TO  REPEAT  THE  PROCESS  ">  ’ 

000 

4ia 

FORMAT  (  /,  5X,  ’  DO  YOU  WANT  TO  MOKE  GRAPH  '•  ’ 

,  X) 

003 

419 

FORMAT  (/,  5X,  ’  DO  YOU  WANT  TO  CORRECT  ->  ’,\) 

004 

450 

FORMAT  (/,  sx,  ’  DO  YOU  WONT  CONTOUR  MOP  ’,\) 

133 


m 


D  Line# 

I  7 

ZLEV 

REAL 

11382 

ZLCW 

REAL 

1 8070 

ZMAX 

REAL 

ians 

ZMIN 

REAL 

18112 

Name 

Type 

Side 

*liCr.;soft  rORTRAN77  V2.  £0  02/8* 


MAIN 

MESHS 

NUMBER 

P3D2D 

PLOT 

PLOTS 

SYMBOL 

WINDOW 

ZCNTUR 

ZLEVEL 


PROGRAM 

SUBROUTINE 

SUBROUTINE 

SUBROUTINE 

SUBROUTINE 

SUBROUTINE 

SUBROUTINE 

SUBROUTINE 

SUBROUTINE 

SUBROUTINE 


1 


Pass  One 


No  Errors  Detected 
£03  Source  Lines 


*  ■  -  “-k  i*-k>  C’Sj.N.lVj.Y.VVj .Mi 


f>j  r.j  nj  PI  nj  nj  nj 


APPENDIX  C 


0  Lino*  l  7 

1  ^STORAGE:  2 

2  »PAGESIZE:38 


Paso 
03-26- 
20 : AG : 

Microsoft  -0RTRAN77  V3.  20  02/ 


THE  PURPOSE  OF  THIS  PROGRAM  IS  TO  COMPUTE  AND  GRAPH  THE  * 
EQUATIONS  OF  ROBERT  P. R0ES3ER  IN  THE  "DISCRETE  STATE-SPACE  * 
MODEL  FOR  LINEAR  IMAGE  PROCESSING".  IT  TRANSFORMS  ALSO  THE  * 
OUTPUT  MATRIX  Y  ACCORDING  TO  FOURIER  ANALYSIS.  * 

• 

evangelos  theofilou  * 

******  *4  **************************************  **********  ******** 


PROGRAM  2D-DATA-FIELD 


*****  VARIABLE  DECLARATIONS  ***** 

REAL  R1  (26,  26)  ,  R2<£6,  26) , SI (26, 26) , S2 (26, 26) , 

*  FR1 (2) , FR2 (2> , TRM ( 4, 4), IV(4),OV(4>, IMGPARI 

CHARACTER*!  ANSWER 


»***  VARIABLE  DECLARATIONS  FOR  PLOTSS  ***** 

CHARACTER*20  CTEXT 

COMMON  /WORK  /Z (26, 26) , ZF(26, 26) , ZL£V(2S) , LDIG(2S> , 

*  LWGT <  26) , MASK ( 3000 ) , VERTEX ( 1 6 ) 

DATA  XLOL/O.  0/,  YLOL/O. 0/, XUPR/S. 3/, YUPR/7. 0/, 

•  ZLOW/1. 0E33/, IPROJ/O/, NRNG/100/ 

*************  MAIN  PROGRAM  **♦»»♦***»**♦ 
*****  ASK  THE  REQUIRED  VALUES  FOR  THE  MODEL  ***** 


*1  * 

10 

WRITE  (*,403) 

wC 

READ  <*, *)  KK 

34 

IF  ( (KK  . LT.  3)  . C 

03 

DO  100  I  ■  1 ,  KK* 1 

36 

DO  100  J  »  1 , KK- 

37 

R1  (1,3)  -  0.  0 

oa 

R2  <  I ,  J)  =•  0.0 

33 

SI  (I,  J)  -  0.0 

40 

S2(  I,  J)  *  0.  O 

41 

42 

100 

CONTINUE 

43 

DO  101  I  *  1,4 

44 

DO  101  J  »  1,4 

43 

TRM(I,J>  *  1 

46 

47 

101 

CONTINUE 

48 

DO  102  1*1,4 

43 

IV  (I)  *  0.0 

30 

OV  ( I )  *  0.  0 

31 

102 

CONTINUE 

C3  U  Oi 


□  Line*  1 
32 
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WRITE  (*,211)  ’  ENTER  INITIAL  CONDITIONS  F OR  HORIZONTAL  Rl't*.  *> 
DO  103  I  =  1,KK 

WRITE  (*,404)  ’Rl(l, ’,!,•):  ’ 


1 

36 

READ  <*,  *)  R1  ( 1,  I ) 

1 

37 

103 

CONTINUE 

38 

39 

WRITE  (*,211)  ’ENTER  INITIAL  COND: 

60 

DO  104  I  »  1 ,  KK 

1 

61 

WRITE  (*,404)  ’ R2  d,’,I,’>:  ’ 

1 

62 

READ  <*,  *)  R2  (1,1) 

1 

63 

104 

CONTINUE 

64 

63 

WRITE  (*,211)  ’ENTER  INITIAL  COND 

66 

DO  103  I  =»  1 ,  KK 

1 

67 

WRITE  (*,403)  ’ SI  (’ ,  I,  ’  ,  1)  :  ’ 

1 

68 

READ  (*, *)  SI ( I,  1 ) 

1 

69 

103 

CONTINUE 

70 

71 

WRITE  (*,211)  ’ENTER  INITIAL  COND 

72 

DO  106  I  =  1,  KK 

1 

73 

WRITE  (*,405)  ’ S£<’ ,  I,  •  ,  1)  :  ’ 

1 

74 

READ  (*, *)  S2 (1,1) 

1 

73 

106 

CONTINUE 

76 

77 

WRITE  (*,211)  ’ENTER  VALUES  FOR  T! 

78 

00(1)  =»  1 

79 

WRITE  (*,409)  ’ bOl s  ’ 

80 

READ  (*, *)  00(3) 

ai 

WRITE  (*,409)  ’ aO 1 :  ’ 

8£ 

READ  <*, *)  0V(4) 

83 

84 

WRITE  (*,211)  'ENTER  ELEMENTS  OF  ' 

83 

TRW (1,2)  -  1 

86 

TRIM  ( 4,  1)  =  1 

87 

WRITE  (*,409)  ’alO:  ’ 

88 

READ  (*, *)  TRM (1,1) 

89 

WRITE  (*,409)  ’ a20 :  ’ 

90 

READ  (*, *)  TRM (2, 1 ) 

91 

WRITE  (*,409)  ’ bl 1 s  ’ 

92 

READ  (*, *)  TEMP 

93 

TRM (1,3)  =  TEMP  +  OV (3) *TRM ( 1 , 1 ) 

94 

WRITE  (*,409)  ’all:  ’ 

95 

READ  (*, *)  TEMP 

96 

TRM (1,4)  a  TEMP  +  00 ( 4 ) *TRM ( 1 , 1 > 

97 

WRITE  (*,409)  ’ b21 :  ’ 

98 

READ  (*, *)  TEMP 

99 

TRM  (2,3)  -  TEMP  +  00  (  3)  *TRM  (2,  1  ) 

100 

WRITE  (*,409)  ’ a2 1 :  ’ 

101 

READ  <*, *)  TEMP 

102 

TRM  (2,  4)  -  TEMP  *  00  ( 4)  *TRM  <2,  1 ) 

■,v,v  v.v.v;.' 


J 


Linen  1 

103 

104 

105 

106 

107 

108 

109 

110 
111 
112 

113 

114 

115 

116 
117 

ns 

119 

120 
121 
122 

123 

124 

125 

126 

127 

128 

129 

130 

131 

132 

133 

134 

135 

136 

137 

138 

139  C 

140 

141 

142 

143 

144 

145 

146 

147 

148 

149 

150 

151 

152 

1  — Jw 


TRM  (4,  3)  =  OV  (3) 
TRM  <4,  4)  *  Qv (4) 
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: 

09-26-8': 
I'J :  42:  3: 
•  20  02/  3* 


WRITE  (*,211)  ’ENTER  VALUES  FOR  THE  INPUT  VECTOR  ( tt  ft)' 

IV  (3)  =*  1 

WRITE  (*,409)  ’ bOO 1  ’ 

READ  (*,*)  I V  < 4 ) 

WRITE  (*,409)  ' blOs  ’ 

READ  <*, *)  TEMP 

IV(1)  »  TEMP  *  IV (4) *TRM ( 1, 1 ) 

IV(2)  =  IV  (4)  *TRM  (2,  1) 

U  »  l .  0 

DO  107  I  =>  1 ,  KK 
DO  107  J  =  1,  KK 

Rl<I  +  l,J>  *  TRM(l,  1)*R1  (I, J>  *  R2(I,3>  +  TRM  < 1 , 3) *S1  ( I , J)  * 

00,  -  ,  ,,  TRM (1,4) *S2 ( I, J)  *  :v(i)*U 

R._(I  +  1, 3)  *  TRM  <2,  1  )  *R1  (  I ,  J)  +  TRm (2.  3) *S1  ( I, 3i  * 

*  TRM <2, 4) *S2 ( I, J)  +  IV(2)*u 

SI  ( I,  3+1 )  =x  (J 

32(1^1)  =*  R1(I,J)  +  OV  (3)  *S1  ( I,  J)  +  OV  (4)  *S2  ( I,  J)  +  IV(4)*U 

107  CONTINUE 

WRITE  (*,205)  ’***♦*  INPUT  VECTOR  ****** 

WRITE  (*,300)  ( IV ( I ) , I  =  1,4) 

WRITE  (*,203)  ’♦****  OUTPUT  VECTOR  ****** 

WRITE  (*,300)  ( OV ( I ) , I  3  1,4) 

WRITE  (*,205)  ’**+**  TRANSITION  MATRIX  ****♦> 

DO  108  I  *  1,4 

WRITE  (*,300)  (TRM ( I , J ) , 3  =  1,4) 

WRITE  (*,210) 

108  CONTINUE 

****  FILL  O’  3  THE  TWO  D I  MENTION  All  GRID  OF  CONTROL  POINTS 
DO  109  I  =  1,26 
DO  109  3  *  l,26 
2(1,3)  *0.0 

109  CONTINUE 

DO  110  I  *  1 , KK 
DO  110  3  *  1 , KK 
Z  (I,  3)  *  R1 ( I,  3) 
no  continue 


WRITE  (*,205) 

DO  111  I  *  1 ,  KK 

WRITE  (*,300)  (Rl(I, 3),  3 
WRITE  (*,210) 


OV (3) *S1  ( I,  3)  +  OV (4) *S2 ( I,  3) 
*****  R 1  MATRIX  ’  ,  KK,  ’  X  ’  ,  KK,  ’ 


1,  KK) 


139 


0  L inert 
1  154 
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continue 


WRITE  ( *,  EOS)  ******  RE  MATRIX  »  ,  KK,  *  X  '  ,  KK,  ’  ****** 
DO  HE  I  =  1,  KK 

WRITE  (*,300)  <R£(I,J>,  J  =  1,KK> 

WRITE  (*,£10) 

CONTINUE 

WRITE  (*,£05)  ******  Si  MATRIX  ’ , KK, ’  X  ’ , KK, ’  ***♦*’ 
DO  113  I  =  1 ,  KK 

WRITE  (*,300)  (S1(I,J),  J  »  1,  KK) 

WRITE  (*,£10) 

CONTINUE 


WRITE  (*, £05)  ’*****  3£  MATRIX  ’  , KK, ’  X  1  , KK, ’  ****** 

DO  114  I  *  1,KK 

WRITE  (*,300)  (S£ ( I, J) ,  J  -  1, KK) 

WRITE  (*,£10) 


172 

114 

CONTINUE 

173 

174 

C 

OUTPUT  THE  Y  MATRIX  ***** 

173 

WRITE 

<  *, 205)  ’ *****  Z 

M  A 

T  R  I  X  ’  ,  KK,  ’ 

176 

WRITE 

<*,  212) 

177 

DO  115 

I  =■  1 ,  KK 

173 

WRIT 

E  (*,300)  (Z(I,J), 

J  = 

1,  KK) 

173 

WRITE  (*,210) 

130 

115 

CONTINUE 

131 

WRITE 

<*, 213) 

132 

133 

WRITE  <*, 413) 

184 

READ 

(*,200)  ANSWER 

135 

IF  <  (ANSWER  .  NE.  ’  Y*  )  . 

AND. 

(ANSWER  .NE.  »y 

136 

137 

c 

*•***♦ 

ASK  THE  PARAMETERS 

FOR 

THE  SRQPH 

133 

20 

WRITE 

(*,  £10) 

133 

WRITE 

<*,  *)  ’ ****  ENT 

E  R 

plot  par 

130 

WRITE 

<*, 410) 

191 

READ 

(*, *)  AZIM 

132 

WRITE 

(*, 411) 

193 

READ 

(*, *)  ELEV 

194 

WRITE 

(*, 413) 

135 

READ 

(*, »)  ITRIM 

136 

WRITE 

(*,414) 

137 

READ 

(*,  *)  ID IV 

138 

WRITE 

(*,  415) 

139 

READ 

(*,  139)  CTEXT 

200 

WRITE 

(*,451) 

£01 

READ 

(*,200)  ANSWER 

£02 

£03 

c 

♦  **■*■» 

INITIALIZE  PL0TS8 

204 

IF  ( (ANSWER  . EO.  ’ Y* )  . 

OR. 

(ANSWER  .  EQ.  ’y’ 

D  une# 
£05 
£06 
£07 
£03 
£0? 
£10 
£11 
£1£ 
£13 
£14 
£15 
£16 
£17 
£18 
£13 
££0 
£21 
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CALL  PLOTS <0,0, £) 

ELSE 

CALL  PLOTS <0, 33, 33) 

END  IF 

CALL  WINDOW ( XLOL, YLOL, XUPR,  YUPR) 

*****  DRAW  THE  MESH  SURFACE  OF  THE  GRAPH  ***** 

CALL  MESHS  <  Z, £6, £6, KK, KK, AZIM, ELEV, 0. 5, 0. 5, a. £5, 6. 5,  IDIV, 0, 

*  3, IPRQJ, 1, ZLOW, 3, ITRIM, MASK, VERTEX) 

*****  ANNOTATION  OF  THE  GRAPH  ***** 

CALL  SYMBOL (1.0, 6. 75, 0. 25, CTEXT, 0. 0, £0) 

CALL  SYMBOL <6. 0, 6.  5, 0. £, ’ £-D  DATA  F I ELD’ , 0. 0,  14 ) 

CALL  SYMB0L<5. 5, 0. 3, 0. £,' AZIMUTHS  ’,0.0,10) 

CALL  NUMBER <333. 0,  333. 0, 0. £, AZIM,  0. 0,  £) 

CALL  SYMBOL  <5. 5, 0. O, O. £, ’ ELEVATION: ’ , 0.  0,  10) 

CALL  NUMBER <333. O, 333. O, O.  £,  ELEV, O.  O,  2) 

DY  =  <Z <  1,  1 ) /30. 0)  *  ELEV 

CALL  P3D2D  < 1 . 0,  1 . 0, Z  < 1 ,  1) -DY, XR, YR> 

CALL  SYMBOL  <  XR, YR,  0. 25, ’ *’ , 0.  O,  1 ) 

CALL  SYMBOL < 1. 0, 0.  1, 0.  £,’ *  =  ORIGIN’ , 0. 0,  10) 

*****  OUTPUT  THE  GRAPH  ***** 

CALL  PLOT <0. 0, 0. 0,  333) 

WRITE  <*, 416) 

READ  <*, £00)  ANSWER 

IF  ( < ANSWER  ,EQ.  ’ Y’  )  .OR.  (ANSWER  . EQ.  ’ y’  > )  GOTO  £0 


£34  21 

WRITE  < *, 418) 

235 

READ  <  *, 200 )  ANSWER 

236 

IF  < (ANSWER  .EQ.  ’ Y’  )  .OR.  (ANSWER 

.  EQ.  ’y’>>  THEN 

237 

£38  C 

****  FILL  0’  3  THE  TWO  DIMENTIONAL 

GRID  OF  CONTROL 

233 

DO  116  I  =  1,26 

1 

£40 

DO  116  J  -  1 , £6 

2 

241 

ZF ( I, J)  *0.0 

2 

£42  116 

CONTINUE 

243 

244 

ZFMAX  »  -3. 3E20 

£45 

ZFMIN  =*  3.  3E20 

246 

DK  *  (KK  -  1)  /  2.  0 

247 

P  =■  3.  141532 

£48 

DO  117  M  »  l ,  KK 

1 

243 

DO  117  N  *  1 , KK 

2 

£50 

rlpart  *  0. 0 

2 

251 

I MG PART  »  0.0 

2 

252 

DO  118  L  =*  1 ,  KK 

3 

£53 

DO  118  K  =  1 , KK 

4 

254 

FR  1(1)  »  COS (-£*P* <L-1 ) * (M 

— DK— 1 ) /KK) 

4 

£55 

FR1  (2)  »  S IN  <  — £*P*  <  L— 1 ) * ( M 

— DK— 1 ) /KK) 

D  Lina#  1 
4  256 

4  257 

4  253 

4  253 

4  260 

4  261 

4  262 

2  263 

2  264 

2  265 

2  266 
2  267 

2  263 

2  263 

2  270 

271 

272  C 

273 

274 

275 
1  276 

1  277 

1  273 


-00 

301  C 

302 

303 

304 

305 

306 
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F R2<1)  =  C0S(-2*P*(K-1)*(N-DK-1) /KK> 

FR2(2>  =  SIN(-2*P*(K-1)*(N-DK-1)/KK> 

RLPART  =  RLPART  +  Z ( L, K > * < FR 1 < 1 ) *FR2 ( l > 

— FR1 (2) *FR2 (2) ) 

IMGPART  =  IMGPART  *  Z  (L,  K>  *  (FR1  ( 1  >  *FR£  (2) 

+FR1 (2) *FR2 <  1 )  > 

CCNTINUE 

ZF(M, N)  »  SORT  (RLPART**2  +  IMGPART**2> 

IF  (ZF(M, N)  .GT.  ZFMAX)  THEN 
ZFMAX  »  ZF(M, N) 

END  IF 

IF  (ZF(M,N>  . LT.  ZFMIN)  THEN 
ZFMIN  »  ZF(M, N> 

END  IF 
CCNTINUE 

*****  OUTPUT  THE  ZF  MATRIX  ***** 

WRITE  <*, 205)  ’ ***  FOURIER  TRANSFORMATION  ' , KK, ’  X  ’ , KK, ’ 
WRITE  <  *,  212) 

DO  113  I  *  1 , KK 

WRITE  (*,300)  <ZF< I,  J) ,  J  =  1 , KK ) 

WRITE  (*,210) 

CONTINUE 


273 

230 

WRITE 

(*, 213) 

281 

WRITE (*, 413) 

2  32 

READ 

<  *, 300) 

ANSWER 

283 

234 

IF  ((ANSWER 

.  NE.  ’  Y’  > 

285  C 

***** 

ASK 

PARAM£~E 

236  30 

WRITE 

(*,  210) 

237 

WRITE 

(*,  *)  • 

***  E  N 

288 

WRITE 

l  *,  410) 

wC*  w 

READ 

(*,  *) 

az  :m 

230 

WRITE 

(*,4U) 

231 

READ 

(*,  *) 

ELEV 

232 

WRITE 

(*, 413) 

233 

READ 

(*,  *) 

ITRIM 

234 

WRITE 

(*,414) 

235 

READ 

(  *,  *) 

IDIV 

236 

WRITE 

(*, 415) 

237 

READ 

(*,  133) 

CTEXT 

233 

WRITE 

(*, 451 ) 

233 

READ 

(*, 200) 

ANSWER 

ER  PLOT  PARAMETERS 


*****  INITIALIZE  PLOT  S3  ***** 


IF  ( (ANSWER  . EQ.  ’ Y’ ) 
CALL  PLOTS (O, O, 2) 
ELSE 

CALL  PLOTS (O,  33,  33) 
END  IF 


(ANSWER  . EQ.  ’  y’  >  >  THEN 


142 


..)3-£2,-g  2 
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WRITE  <*,  4£0) 

READ  (*,£00)  ANSWER 


CALL  WINDOW (XLOL, YLOL, XUPR, YUPR) 


IF  < (ANSWER  .ED.  • Y’  )  .OR.  (ANSWER  . EQ.  ’y’>>  THEN 
DLEV  =  (ZFMAX-ZFMIN) /FLOAT (KK) 

CALL  ZLEVEL  ( ZF,  £6,  £6,  KK,  KK,  DLEV,  ZLEV,KK*1) 

DO  136  I  =  1 , KK+ 1 
LDIO(I)  -  £ 

LWGT ( I )  -  1 
CONTINUE 

CALL  ZCNTUR (ZF,  £6,  £6,  KK,  KK,  0.  5, O. 5, 3. £5, 6. 5, ZLEV, LDIG, LwGT, 
KK+1 , 0. 10, 10) 

CALL  SYMBOL  <3.  5,0.  0,0.  £,  ’CONTOUR  MAP’,  0.0,  11) 

ELSE 

*****  DRAW  THE  MESH  SURFACE  OF  THE  GRAPH  ***** 

CALL  MESHS  <  ZF , £6, £6, KK, KK, AZIM,  ELEV, 0.  3,  0. 5,  3. £3,  6.  3,  IDIV, 0, 
3, I PRO J, 1, ZLOW, 3, ITRIM, MASK, VERTEX) 


*****  ANNOTATION  OF  THE  GRAPH  ***** 

CALL  SYMBOL (3. 3,  0. 3,  0.  £,’ AZIMUTH:  ’,0.0,10) 

CALL  NUMBER  <399. 0, 993. 0, 0.  £,  AZIM, 0. 0,  £> 

CALL  SYMBOL (3. 3, 0. 0, 0.  £,  ’ ELEVATION: • , 0. 0,  10) 
CALL  NUMBER <999.  0, 399.  0,  O.  £,  ELEV,  0.  0,  £) 

DY  =  (ZF<1, l)/90. 0)  *  ELEV 

CALL  P3D£D  <1.0,  1.0,  ZF  ( 1 ,  1 ) -DY,  XR,  YR) 

CALL  SYMBOL  <  XR,  YR, 0. £5,  ’ *’  , 0. 0,  l ) 

CALL  SYMBOL ( 1. 0, 0.  1, 0. £,  ’ *  »  ORIGIN’ , 0. 0,  10) 
END  IF 

CALL  SYMBOL (1.0, 6.  73,  0. £3, CTEXT, 0. 0,  £0) 

CALL  SYMBOL (6. 0, 6. 3, 0. £, ’ £— D  DFT’,0.0, 7) 


*****  OUTPUT  THE  GRAPH  ***** 

CALL  PLOT  <0.  0, 0. 0,  999) 

WRITE  (*,416) 

READ  <*, £00)  ANSWER 

IF  ((ANSWER  . £□.  ’ Y’  )  .OR.  (ANSWER  .  EQ.  ’ y’  ) )  GOTO  30 

END  IF 

WRITE  (*,417) 

READ  (*,£00)  ANSWER 

IF  < (ANSWER  .EQ.  ’ Y’ >  .OR.  (ANSWER  .ED.  ’ y’ ) )  GOTO  10 
STOP 


199  FORMAT <A£0) 

£00  FORMAT (A) 

£03  FORMAT ( / ,  18X,A£3,  I£, A3,  I£, AS,  /> 

£10  FORMAT ( ) 

£11  FORMAT </,3X,A36) 

£1£  FORMAT(/,£X,  ’  (AZIMUTH  3EO. 0) ’ , 46X, ’  (AZIMUTH  £30.0)’,/) 


D  Line* 
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353 

213 

FORMAT!/,  2X,  ’ 

(AZIMUTH  050.  0) ’, 46X,  ’  (AZIMUTH  140.0)’,/) 

353 

300 

FORMAT ( 10 (F7. 

2,  IX)  ) 

360 

400 

FORMAT  OX,  \  ) 

361 

403 

FORMAT ( / , SX, ’ 

DIMENSION  OF  OUTPUT <K-lto20> .  ’ , \ > 

363 

404 

FORMAT  <  5X , AS, 

12,  A3,  \> 

363 

405 

FORMAT  (SX,  A3, 

12,  AS,  \) 

364 

406 

FORMAT  <SX, A2, 

12,  A4,  \> 

365 

407 

FORMAT (5X, A2, 

12,  12,  A3,  \> 

366 

408 

FORMAT  <SX,  A3, 

12,  A3,  \) 

367 

403 

FORMAT (5X,  AS, 

\) 

363 

410 

FORMAT (/, SX ,  » 

AZ IMUTH ( 0. 0  to  360.0  DEGREES):  ’,\> 

363 

41 1 

FORMAT  </,  5X,  ’ 

ELEVATION (30. 0  to  -30.0  DEGREES): 

370 

413 

FORMAT (/,  SX,  ’ 

TRIM(0«N0,  l*Xs, 2»Ys)  :  ’,\) 

371 

414 

FORMAT ( / , SX, • 

2, 4  OR  a  SUBGRIDS:  ’ , \) 

372 

415 

FORMAT ( / , SX,  ’ 

TITLE  OF  GRAPH (UP  TO  20  CHAR) :  ’,\) 

373 

416 

FORMAT  </, 5X , ’ 

DO  YOU  WANT  TO  CHANGE  PARAMETERS'  ’,\) 

374 

417 

FORMAT ( / ,  SX,  ’ 

DO  YOU  WANT  TO  REPEAT  THE  PROCESS'  ’,\) 

375 

413 

FORMAT ( / , SX, » 

DO  YOU  WANT  FOURIER  TRANSFORMATION  t>  .  (N) 

376 

413 

FORMAT!/, SX, ’ 

DO  YOU  WANT  TO  MAKE  GRAPH  '  ’,\) 

377 

420 

FORMAT!/,  SX,  ’ 

DO  YOU  WANT  CONTOUR  MAP  7  >t\) 

373 

451 

FORMAT !/, SX, ’ 

SEND  GRAPH  TO  THE  PRINTER ( Y  or  N) :  >,\) 

373 

END 

Name 

Type 

Offset  P 

Class 

ANSWER 

CHAR* 1 

1 1060 

AZIM 

REAL 

11062 

COS 

INTRINSIC 

CTEXT 

CHAR*20 

11074 

DK 

real 

11114 

DLEV 

REAL 

11163 

DY 

REAL 

1 1034 

ELEV 

REAL 

1 1 066 

FLOAT 

INTRINSIC 

FRl 

REAL 

10882 

r  R«£ 

REAL 

10830" 

i 

INTEGER*2 

10356 

IDIV 

INTEGER*2 

11072 

imgpar 

REAL 

11142 

I  PRO  J 

INTEGER*2 

10350 

ITRIM 

INTEGER*2 

11070 

IV 

REAL 

10338 

J 

INTEGER*2 

10364 

K 

INTEGER*2 

1 1 154 

KK 

INTEGER*2 

10354 

L 

INTEGER*2 

11146 

LDIG 

INTEGER*2 

5512 

/WORK  / 

LWGT 

INTEGER*2 

SS64 

/WORK  / 

M 

INTEGER*2 

1 1 122 

MASK 

INTEGER*2 

5616 

/WORK  / 

N 

INTEGER*2 

11130 

,  %  *  *• 

m 

■ss 


MS 

V  V 

» v  » 

M 


».  » >  » 


144 


*.-<«*••  o  •\_’<  «  * 


,'  / 

■  ■ .  ' 


L  lV/- 


3 .  '■ 


M  *  cr 


D  Lina#  1  7 


NRNG 

INTEGER*^ 

1 0352 

OV 

REAL 

10314 

P 

REAL 

11118 

R1 

REAL 

R£ 

REAL 

2706 

RLPART 

REAL 

11138 

SI 

REAL 

3410 

S£ 

REAL 

8114 

SIN 

INTRINSIC 

SORT 

INTRINSIC 

TEMP 

REAL 

10396 

trm 

REAL 

10818 

U 

REAL 

1 1 000 

VERTEX 

REAL 

11616 

/WORK  / 

XLOL 

REAL 

10330 

XR 

REAL 

1 1038 

XUPR 

REAL 

10338 

YLOL 

REAL 

10334 

YR 

REAL 

11102 

YUPR 

REAL 

10342 

Z 

REAL 

0 

/WORK  / 

ZF 

real 

2704 

/WORK  / 

ZFMAX 

REAL 
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ZFMIN 

REAL 
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the  purpose  of  this  program  is  to  cgmpute  and  graph  -he 
EQUATIONS  OF  ROBERT  P. RQESSER  IN  THE  "DISCRETE  STATE— SPACE 
MODEL  FOR  LINEAR  IMAGE  PROCESSING".  IT  TRANSFORMS  ALSO  THE 
OUTPUT  MATRIX  Y  ACCORDING  TO  FOURIER  ANALYSIS. 

EVANGELOS  THEOFILOU 

I 

PROGRAM  2D— DAT A -FI ELD 

*****  VARIABLE  DECLARATIONS  ***** 

REf,L  *  <26,  £S,  4)  ,  31  (£S,  28,  4)  ,  S2  (££,  2S  .  4)  , 

*  R1  (2)  ,  R£  <£) ,  TRM  <12,  12),  IV(12),QV<1£>,  IMGPART 
CHARACTER*!  ANSWER 

****  VARIABLE  DECLARATIONS  FOR  PLQTaa  ***** 

CHARACTER *20  CTEXT 

COMMON  /WORK  /Z  <26, 26) , ZF(2S, 26) , ZLEV<£S> , LDIG (£6) , 

*  LWGT ( 26 ) , MASK (2000)  ,  VERTEX  <1S) 


DATA 


XLOL/O. 0/, YLCL/O.  0/, XUPR/8.  3/,  YUPR/7.  O/, 
ZLQW/ 1. 0E23/, IPROJ/O/, NRNG/ 100/ 


*************  MAIN 


P  R  0  G  R  A 


*****  aSK  THE  REQUIRED  VALUES  FOR  THE  MODEL  ***** 
10  WRITE  <*, 401 ) 

READ  <*, *)  N 

IF  < <N  .LT.  1)  .OR.  ( N  .GT.  4))  GOTO  10 

2  WRITE  <  *,  402) 

READ  <*, *)  M 

IF  ( <M  .LT.  1)  .OR.  <M  . GT.  4))  GOTO  2 

3  WRITE  <*,  403) 

READ  <*,  *)  KK 

IF  <KK  .GT.  23)  GOTO  3 

DO  100  I  »  1 , KK* 1 
DO  100  J  =  1 , KK* l 
DO  1 00  L  *  1 , N 
R  <  I,  J,  L)  =  0.0 

51  <  I,  J,  L)  =*  O.  0 

52  ( I,  J,  L)  =»  0.0 
100  CONTINUE 

DO  101  I  =  1,N*2»M 
DO  101  J  =  1,N+£*M 
TRM  <  I  ,  J)  »  0.0 


l 


Ul  -\J  OJ 


0  Line**  1 


CONTINUE 

DC  102  I  =  1,N+2*M 
IV (I)  =  0.0 
OV  <  I )  =0.  O 
CONTINUE 
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WRITE  (*,211)  'ENTER  INITIAL  CONDITIONS  FOR  HORIZONTAL  R(».#)’ 
DO  103  I  =  1,KK 
DO  103  J  =■  1 ,  N 

WRITE  (*,404)  ’  R’  ,  J,  ’  (1,  ’  ,  I,  ’  )  ;  ’ 

READ  (*, *)  R( 1,  I,  J> 

CONTINUE 

WRITE  (*,211)  'ENTER  INITIAL  CONDITIONS  FOR  VERTICAL  SI (#.*»>  ’ 
DO  104  I  =  1 , KK 
DO  104  J  =  1,M 

WRITE  (*,405)  ’  31  <’  ,  J,  ’  )  (’  ,  I,  '  ,  1 )  :  ’ 

READ  (*,  *)  Si  <  I,  1,  J) 

CONTINUE 

WRITE  (*,211)  'ENTER  INITIAL  CONDITIONS  FOR  VERTICAL  S2(#. #>  ’ 
DO  105  I  =*  1 ,  KK 
DO  105  J  =  1 , M 

WRITE  (*,405)  ’  S2C  ,  J,  ’  )  <’  ,  I,  ’  ,  1)  !  ’ 

READ  (*, *)  S2( I,  1,  J) 

CONTINUE 

WRITE  (*,211)  'ENTER  VALUES  FOR  THE  INPUT  VECTOR(*.  *)  ’ 

IV  (1)  *  1.0 
DO  1 06  I  *  1 ,  M 

WRITE  (*,408)  ’ a ( 0 • , I , ’ ) :  ’ 

READ  (*, *)  IV(N+I) 

TRM (N+I,  N+l )  =  —IV (N+I ) 

CONTINUE 
DO  107  I  =  1,01 

WRITE  (*,408)  ’ b (O’  ,  I ,  ’  )  s  ’ 

READ  (*, *)  IV(N+M+t) 

TRM(N+M+I,  N+l)  =  -IV(N+M+I) 

CONTINUE 

DO  108  I  =  1,M-1 

TRM(N+I,N+I+1)  =  1.0 
CONTINUE 
DO  109  I  =»  1 ,  M—  1 

TRM(N+M+I, N+M+I+l )  =  1.0 
CONTINUE 


WRITE  (*,211)  'ENTER  ELEMENTS  OF  THE  TRANSITION  MATRIX!#.#)  ’ 
DO  110  I  =*  1 ,  N 

WRITE  (*,406)  ' a ( ’ ,  I , ’ O  >  :  ’ 

READ  (*, *)  TEMP 
TRM ( 1 , I )  =  -TEMP 
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no  Continue 

TRM  (  1 ,  N+- 1  )  »  -1.0 
DO  111  I  =  2,  N 

TRW  (1,1-1)  =  1.0 

111  CONTINUE 

DO  112  I  =  i,n 
DO  112  J  =  1 , N 

WRITE  (*,407)  1  a(’,J,  I,  '  )  :  ’ 

READ  (*,  *>  TEMPI 

TRM  (  I +N,  J )  a  TEMPI  *  TRM(1,J)  *  IV(N+I> 

112  CONTINUE 

DO  113  I  =  l.M 
DO  113  J  =  1,N 

WRITE  (*,407)  ’  b  ( ’  ,  J,  I,  *  )  :  ’ 

READ  (*, *)  TEMPI 

TRM ( I+N+M, J)  =«  TEMPI  *  TRM  (  1 ,  J )  *  IV(N+M-*-I) 

113  CONTINUE 

WRITE  (*,211)  ’ENTER  VALUES  FOR  THE  OUTPUT  VECTOR <#.#) 

WRITE  (*,409)  ’b(OO):  ’ 

READ  (*, *>  TEMP 
QV (N+l )  =  -TEMP 
OV  ( N+M+  1  )  =  1.0 
DO  114  I  =*  l,  N 

WRITE  (*,406)  *  b  ( *  ,  I ,  ’  0)  :  ’ 

READ  (*, *)  TEMPI 

QV(I)  =  TEMPI  *  TRM (1,1)  *  TEMP 

114  CONTINUE 

U  =  1.0 

DO  115  I  =»  1 ,  KK 
DO  115  J  =  1 , KK 

DO  US  II  =  l,N*2*M 

IF  (II  . LE.  N)  THEN 
DO  117  JJ  =  1 ,  N-*-2*M 

IF  (JJ  . LE.  N)  THEN 

R(I*1,  J,  II)»R(I  +  1,  J,  I  I )  -+-TRM  (  I  I ,  J  J )  *R  (  I ,  J,  J  J ) 

END  IF 

IF  ((JJ  . GT.  N)  .AND.  (JJ  . LE.  N+M) )  THEN 

R(I*1,  J,  II)=R(I*1,  J,  II)*TRM(II,JJ)*S1  (I,  J,  JJ-N) 
END  I F 

117  CONTINUE 

R  (  1*1,  J,  I  I ) -R  (  1*1,  J,  1 1  )  *  IV(II)  *  U 
END  IF 

IF  ((II  . GT.  N)  .AND.  (II  . LE.  N+M) )  THEN 
DO  llfl  JJ  =  1,N+2*M 

IF  (JJ  . LE.  N)  THEN 

SI  (  I,  J*l,  I I-N)  =»  SI  (  I,  J+l,  I  I-N)  *  TRM  (  I  I ,  JJ )  * 

*  R ( I , J, JJ ) 
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END  IF 

IF  ((JJ  .  3T.  N).AND.  <JJ  .  LE.  N+M)  )  THEN 

SI  (I,  J+l,  II-N)  =  SI  (I,  J+l,  II-IM)  +  TRM(II,JJ>* 

*  S1CI,J,  JJ-N) 

END  IF 

118  CONTINUE 

SI  (I,  J+l,  II-N)  -  SI  (I,  J+l,  II-N)  +  IV(II)  *  U 
END  IF 

IF  (II  .  GT.  N+M)  THEN 
DO  113  JJ  =  1,N+£*M 

IF  <JJ  .LE.  N)  THEN 

S£(I,  J+l,  II-N-M)  =  S2<I,  J+l,  II-N-*!)  *  TRM(II.JJ) 

*  *  3(1,  J, JJ) 

END  IF 

IF  (<JJ  . GT.  N>  .AND.  (JJ  . LE.  N+M > )  THEN 

S£( I, J+l, II-N-M)  =  S2(I,  J+l,  II-N-M)  +  TRM  <  1 1 ,  J  J ) 

*  *  ’SI  (I,  J,  JJ-N) 

END  IF 

IF  (JJ  . GT.  N+M)  THEN 

S2( I, J+l,  II-N-M)  -  S£(I, J+l, II-N-M)  +  TRM ( 1 1 , JJ ) 

*  *  S£ ( I, J, JJ-N-M) 

END  IF 

113  CONTINUE 

S2 ( I,  J+l,  I I-N-M)  *  S£(I, J+l,  II-N-M)  ♦  I V ( I  I >  *  U 

END  IF 

116  CONTINUE 

U  =  0.  0 
113  CONTINUE 

WRITE  (*,203)  <*+***  INPUT  VECTOR  ♦****> 

WRITE  (*,300)  (IV  (I),  I  =  1,N+£*M> 

WRITE  (*,203)  >****♦  OUTPUT  VECTOR  ****+> 

WRITE  (*,300)  <OV(I),I  =  1,N+2*M) 

WRITE  (*,205)  ’**+**  TRANSITION  MATRIX  ***+*» 

DO  120  I  =  1 ,  N+2*M 

WRITE  (*,300)  (TRM(I,J),J  =  1,N+2*M) 

WRITE  (*,210) 

120  CONTINUE 

****  FILL  O’ s  THE  TWO  DIMENTIONAL  GRID  OF  CONTROL  POINTS  **** 

DO  121  I  =  1,26 
DO  121  J  =  1,26 
Z(I,J)  =0.0 

121  CONTINUE 


DO  122  I  =  1 ,  KK 
DO  122  J  =  1 , KK 

DO  123  LL  =  1 ,  N+2+M 


-  V* ka  U  /* *  ^ 


D 

Line#  1 

7 

3 

305 

IF  (LL  . L£.  N)  ~h£N 

3 

206 

Z  (I,  J)  =  Z  ( I,  J)  +  00  (LL)  * 

3 

307 

END  IF 

3 

£08 

IF  <<LL  . GT.  N).AND. (LL  . LE. 

3 

303 

Z  ( I,  J)  =  Z ( I, J)  -<-00  (LL)  * 

3 

210 

END  I F 

vj 

211 

IF  (LL  .  GT.  N+M>  THEN 

Ml 

213 

Z  (I,  J)  =*  Z  ( I,  J)  +  00  (LL)  * 

3 

213 

END  IF 

j 

214 

123 

CONTINUE 

a 

313 

122 

CONTINUE 

316 

317 

WRITE  (*,205)  ’*****  R  M  A  T  R  I 

218 

DO  124  I  *  1 , KK 

i 

313 

DO  125  L  =  1,N 

2 

220 

WRITE  (*,300)  <  R  ( I ,  J,  L)  ,  J  =*  1 

2 

221 

125 

CONTINUE 

1 

222 

WRITE  (*,210) 

1 

323 

124 

CONTINUE 

224 

225 

WRITE  (*,205)  ’*****  SI  M  A  T  R  I 

226 

DO  126  I  =•  1 ,  KK 

1 

227 

DO  127  L  -  1, M 

o 

228 

WRITE  (*,300)  (S1(I,J, L),  J  =■ 

3 

223 

127 

CONTINUE 

1 

230 

WRITE  (*,210) 

1 

231 

126 

CONTINUE 

232 

333 

WRITE  (*,205)  >***♦*  32  M  A  T  R  I 

234 

DO  128  I  =  1 ,  KK 

1 

235 

DO  123  L  »  l,M 

2 

236 

WRITE  (*,300)  (S3(I,J, L),  J  * 

a 

237 

133 

CONTINUE 

i 

238 

WRITE  (*,210) 

i 

233 

128 

CONTINUE 

240 

241  C 

*****  OUTPUT  THE  Z  MATRIX  ***** 

242 

WRITE  (*,205)  ’*****  Z  M  A  T  R  I 

243 

WRITE  (*,212) 

244 

DO  130  I  *  1 ,  KK 

i 

245 

WRITE  (*,300)  (Z(I,J),  J  =  1 ,  KK ) 

i 

246 

WRITE  (*,210) 

i 

247 

130 

CONTINUE 

248 

WRITE  (*,213) 

243 

250 

WR I TE ( *,  413) 

251 

READ  (*,200)  ANSWER 

Micr-osor'5  FQR7RAN77 


IF  ((ANSWER  . NE.  ’ Y»  )  .AND.  (ANSWER  .  NE.  ’ 


354  C 

355  30 


*****  psk  THE  PARAMETERS  FOR  THE  GRAPH 
WRITE  (*,310) 


y’ ) )  GOTO  21 
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281 
232 

283  C 

284 

285 

286 
237 
288 
283 

230 

231 

232 

233 

234 

235  C 

236 

237 

238 
233 

300 

301  21 

302 

303 

304 

305  C 
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WRITE  (*,  *)  ’  ****  ENTER  PLOT  PARAMETERS  ***•< 

WRITE  (*,410) 

READ  <*, *)  AZIM 
WRITE  (*,411) 

READ  (*, *)  ELEV 
WRITE  (*,413) 

READ  (*, *)  ITRIM 
WRITE  (*,414) 

READ  (*, *)  IDIV 
WRITE  (*,415) 

READ  (*, 133)  CTEXT 
WRITE  (*,451) 

READ  (*,200)  ANSWER 

*****  INITIALIZE  PL0TS8  ***** 

IF  ((ANSWER  . EQ.  ’ Y’ )  .OR.  (ANSWER  . EQ.  ’ y’ ) >  THEN 
CALL  PLOTS (0,0,2) 

ELSE 

CALL  PLOTS (O,  39,  33) 

END  IF 

CALL  WINDOW ( XLOL,  VLOL,  XUPR, YUPR) 

*****  DRAW  THE  MESH  SURFACE  OF  THE  GRAPH  ***** 

CALL  MESHS  <Z,  26.  26,  KK,  KK,  AZIM,  ELEV, 0. 5, O. 5,  8. 25, 6. 5,  IDIV, 0. 

►  3, IPROJ, 1, ZLOW, 3, ITRIM, MASK, VERTEX) 

*****  ANNOTATION  OF  THE  GRAPH  ***** 

CALL  SYMBOL ( 1. 0, 6. 75, 0. 25, CTEXT, 0. 0, 20) 

CALL  SYMBOL (6. 0,  6.  5,  O.  2,  ’ 2-D  DATA  FIELD* , 0. O,  14) 

CALL  SYMBOL (5. 5,  0.  3,  0.  2,  *  AZIMUTH:  ’,0.0,10) 

CALL  NUMBER (333.  0,  333.  0,  0. 2, AZIM, 0.  O, 2) 

CALL  SYMBOL (5. 5.  0.  0,  0.  2,  ’ ELEVATION:’ , 0. 0,  10) 

CALL  NUMBER (333.  0,  393. O, O.  2, ELEV, 0. 0, 2) 

DY  »  (Z(l,  D/30.0)  *  ELEV 

CALL  P3D2D (1.0,  1.0, Z  < 1,  1 ) -DY, XR, YR) 

CALL  SYMBOL (XR, YR, 0.  25, ’ *’ , O. 0,  1 ) 

CALL  SYMBOL  (  1. 0,  0.  1,  O.  2,  ’  *  =*  ORIGIN’  ,  0.  0,  10) 

*****  OUTPUT  THE  GRAPH  ***** 

CALL  PLOT (0. 0, 0. 0, 393) 

WRITE  (*,416) 

READ  (*,200)  ANSWER 

IF  ((ANSWER  .EQ.  ’ Y’ )  .OR.  (ANSWER  . EQ.  ’ y’ > >  GOTO  20 

WRITE (*,  418) 

READ ( *, 200 )  ANSWER 

IF  ((ANSWER  .EQ.  ’ Y’  )  .OR.  (ANSWER  . EQ.  ’y’)>  THEN 

****  FILL  0’ s  THE  TWO  DIMENTIONAL  GRID  OF  CONTROL  POINTS  **** 
DO  132  I  »  1,26 
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DO  135  J  =  1.56 
ZF(I,J)  =0.0 
CONTINUE 

ZFMAX  =  -9.  9E20 
ZFMIN  =  9. 9E50 
DK  =  (KK  -  1)  /  2.  O 

P  =  3. 141592 
DO  133  M  =  1 ,  KK 
DO  133  N  =  1 , KK 
RLPART  =  0.  0 
IMGPflRT  =  0.  O 
DO  134  L  =  1,KK 
DO  134  K  =  1,KK 

R 1  (  1 )  =  COS (-2*0* (L-l > * (M-DK-1 ) /KK) 
R1 <2)  =  SIN (-2*P*(L-1 ) *(M-DK-1 > /KK) 
R2  <  1 )  =  COS (-2*0* (K-l >* (N-DK-1 > /KK) 


R2 (2)  =  SIN  <— £*P*  <K— 1 ) *  <N— DK— 1 ) /KK) 
RLPART  =  RLPART  +■  Z  ( L,  K >  *  (  R 1  ( 1 )  *R, 

— R 1  (5)  *R, 

IMGPflRT  =  IMGPART  *  Z  (L,  K>  *  (  R1  ( 1  >  *R 

+R1  (2)  *R, 

CONTINUE 

ZF  ( M,  N )  =  SORT <RLPART**2  +  IMGPART**2) 
IF  (ZF(M, N)  .GT.  ZFMAX )  THEN 
ZFMAX  =  ZF(M,N) 

END  IF 

IF  <ZF(M, N>  . LT.  ZFMIN)  THEN 
ZFMIN  =  ZF<M,  N> 

END  IF 
CONTINUE 


IMGPflRT 


2(1) 
2(2)  ) 
2(2) 
2(D) 


*****  OUTPUT  THE  ZF  MATRIX  ***** 

WRITE  (*,203)  ’ ***  FOURIER  TRANSFORMATION  ’  ,  KK,  ’  X  ’ , KK, '  ***’ 

WRITE  (*,212) 

DO  133  I  =  1 ,  KK 

WRITE  (*,300)  ( ZF ( I , J ) ,  J  =  1,  KK) 

WRITE  (*,210) 

CONTINUE 
WRITE  (*,213) 

WRITE ( *,  419) 

READ  (*,200)  ANSWER 

IF  ((ANSWER  .NE.  •  Y*  >  .AND.  (ANSWER  . NE.  ’ y' ) )  GOTO  22 


*****  ASK  THE  PARAMETERS  FOR  THE  GRAPH  ***** 
WRITE  (*,210) 


WRITE  (* 
WRITE  <* 
READ  (* 
WRITE  (* 


•)  ’ •**  ENTER  PLOT  PARAMETERS  *»*’ 

410) 

*)  AZIM 
411  ) 


•*  •  .  *  .k*  k  r.m 


(n 
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7 

338 

READ 

(*, *i  EL2V 

353 

WRITE 

(*,  413) 

360 

READ 

(*,  *)  ITRIM 

361 

WHITE 

(*,  414) 

362 

READ 

(*,*)  IDIV 

363 

WRITE 

(*,  413) 

364 

READ 

(*,  199)  CTEXT 

363 

WRITE 

(*,  431) 

366 

367 

READ 

(*,200)  ANSWER 

368  C 

***** 

INITIALIZE  PLOT 88  ***** 

369 

IF  ((ANSWER  ,  EQ.  ’  Y’  )  .OR.  (1 

370 

CALL 

PLOTS (0, 0, 2) 

371 

ELSE 

372 

CALL 

PLOTS (0, 39, 39> 

373 

374 

END  IF 

373 

WRITE 

(*,  420) 

376 

READ  (*,200)  ANSWER 

Pane  d 
09-25-92 
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377 

378 

379 

380 

381 

382 

383 

384 
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386 

387 

388 

389 

390 

391  C 

392 

393 

394 
393  C 

396 

397 

398 

399 

400 

401 

402 

403 

404 

403 

406 

407 

408  C 


CALL  WINDOW < XLOL,  YLOL,  XUPR,  YUPR) 


IF  ((ANSWER  .  EQ.  » Y’  >  .OR.  (ANSWER  .  EQ.  * y’  > )  THEN 
DLEV  *  ( ZFMAX— ZFMIN) /FLOAT (KK) 

CALL  ZLEVEL(ZF,  26,  26,  KK,  KK,  DLEV,  ZLEV,  KK+1) 

DO  136  I  =  l.KK+1 
LOIG(I)  -  2 
LWGT ( I )  =1 

CONTINUE 

CALL  ZCNTUR ( ZF , 26, 26, KK, KK, 0. 3, 0. 3, 8. 23, 6. 5, ZLEV, LDIG, LWGT, 
KK*1,  0.  10,  10) 

CALL  SYMBOL (3. 3,  0.0,  0.2,  ’CONTOUR  MPP’,0. O,  11> 

ELSE 

*****  DRAW  THE  MESH  SURFACE  OF  THE  GRAPH  ***** 

CALL  MESHS ( ZF, 26, 26,  KK, KK, AZIM, ELEV, 0. 3, 0. 3, 8. 25, 6. 5,  IDIV, 0, 
3, IPROJ, 1, ZLOW, 3, ITRIM, MASK, VERTEX) 


*****  ANNOTATION  OF  THE  GRAPH  ***** 

CALL  SYMBOL (3.  5,  0.  3,  0.2,  •  AZIMUTH:  ’,0.0,10) 

CALL  NUMBER (999. 0, 999. 0,  0.  2,  AZIM,  0.  0,  2) 

CALL  SYMBOL (3. 3, 0. 0,  0.  2,  ’  ELEVATION:  ’  ,  0.  0,  10) 
CALL  NUMBER (999.  0,  999.  0,  0.  2,  ELEV, 0.  0,  2) 

DY  »  (ZF ( 1,  1 ) /90. 0)  *  ELEV 

CALL  03D2D < 1. 0,  1. 0, ZF ( 1,  1 )  -DY,  XR,  YR> 

CALL  SYMBOL  <  XR, YR, 0. 23, ’ *’  , 0. 0,  1 ) 

CAL_  SYMBOL ( 1. O, O. 1, O. 2, ’ *  -  OR IGI N’ , 0. 0, 1 0 ) 
END  IF 

CALl  SYMBOL (1.0, 6. 73, 0. 23, CTEXT,  0. O,  20) 

CALL  SYMBOL (6.  0,6.  3,  0.  2,  ’2-0  DFT’, 0.0,7) 


OUTPUT  THE  GRAPH 


03-26-8 
13: 32:0 


D  Lire# 

1 
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403 

CALL  PLOT (0. 0, O. 0, 333) 

410 

WRITE  (*,416) 

41  1 

READ  (*,200)  ANSWER 

412 

IF  ((ANSWER  . EO.  ’  Y’  )  .OR.  (ANSWER  . EQ.  ’y’)>  GOTO  30 

413 

£2 

END  I F 

414 

WRITE  (*,417) 

415 

READ  (*,200)  ANSWER 

416 

IF  ((ANSWER  .EQ.  ’ Y’ )  .OR.  (ANSWER  . EQ.  ’ y’ ) )  GOTO  10 

417 

STOP 

418 

413 

199 

FORMAT (AGO) 

420 

200 

FORMAT (A) 

421 

203 

FORMAT*/,  1SX, A23,  12, A3,  12, AS, /) 

422  £10  FORMAT*) 

423  £11  FORMAT  <  / , 5X,  A46 ) 

424  £  1  £  FORMAT (/, 2X,  ’  (AZIMUTH  320. 0) ’, 46X, ’ (AZIMUTH  £20.0)’,/) 

425  £13  FORMAT*/, £X,  ’  (AZIMUTH  030.  0)  ’  ,  46X,  M AZ IMUTH  140.0)’,/) 

4£6  300  FORMAT ( 10 (F7.  £,  IX)  ) 

4£7  400  FORMAT  OX,  \) 

4£S  451  FORMAT*/, 5X, ’ SEND  GRAPH  TO  THE  PRINTER ( Y  or  N) :  ’,\) 

4£3  4<.n  FORMAT  (/,  5X,  ’  NUMBER  OF  HORIZONTAL  STATES < N= 1 t o4 ) :  ’,\) 

430  402  FORMAT*/, 5X, ’ NUMBER  OF  VERTICAL  STATES ( M= 1 to4 ) :  >,\) 

431  403  FORMAT (/, 3X,  ’  DIMENSION  OF  OUTPUT ( lto£5)  :  ’  ,  \ ) 

432  404  F0RMAT(3X,  Al,  I£, A3,  12, A3,  \> 

433  405  FORMAT <5X,  A3,  I£, A£,  12, A5, \) 

434  406  FORMAT  <5X,  A2,  12,  A4,  \) 

425  407  FORMAT (3X, A2, I£, 12, A3, \) 

436  40S  FORMAT (3X, A3,  12,  A3,  \ ) 

437  403  FORMAT (SX, A8, \ ) 

438  410  FORMAT*/,  3X, ’ AZIMUTH(0. 0  to  360.0  DEGREES):  ’,\> 

433  411  FORMAT*/, SX,  ’ ELEVATION <30. 0  to  -30.0  DEGREES):  ’,\) 

440  413  FORMAT (/, 5X,  ’ TRIM(0=N0,  l  =  Xs, 2»Ys)  :  ’,\) 

441  414  FORMAT  </,5X,  ’£,4  OR  a  SUSG.RIDS:  ’,\) 

442  413  FORMAT*/, 5X, ’ TITLE  OF  GRAPH (UP  TO  20  CHAR):  ’,\) 

443  416  FORMAT*/,  5X,  ’  DO  YOU  WANT  TO  CHANGE  PARAMETERS'1  ’,\> 

444  417  FORMAT*/,  3X,  ’  DO  YOU  WANT  TO  REPEAT  THE  PROCESS?  ’,\> 

445  418  FORMAT*/,  5X,  ’  DO  YOU  WANT  FOURIER  TRANSFORMAT  I  ON  -1  ’,\) 

446  413  FORMAT*/,  5X,  ’  DO  YOU  WANT  TO  MAKE  GRAPH  ?  ' , 

447  420  FORMAT </,  5X,  ’  DO  YOU  WANT  CONTOUR  MAP  ?  ’,\) 

448  END 


Name 

Tyoe 

Offset  P 

C  1  ckS3 

ANSWER 

CHAR* 1 

33434 

AZIM 

REAL 

33426 

COS 

CTEXT 

CHAR*20 

33448 

INTRINSIC 

DK 

REAL 

33488 

DLEV 

REAL 

33336 

DY 

REAL 

33468 

ELEV 

REAL 

33440 

154 


»/  \  •.  \  %  *, 


iX»  (i|  a) 


0  Llnstt 
FLOAT 


I 

INTEGERS 

33163 

IDIV 

INTEGSR*2 

33446 

II 

INTEGER*2 

33336 

IMPPAR 

REAL 

33512 

I  PRO  J 

INTEGER*2 

33153 

ITRIM 

INTEGERS 

33444 

IV 

REAL 

33042 

J 

INTEGER*2 

33176 

JJ 

INTEGER*2 

33344 

K 

INTEGER*2 

33522 

KK 

INTEGER*2 

33 1 66 

L 

INTEGER*2 

33134 

LDIG 

INTEGER*2 

5512 

LL 

INTEGER*2 

33334 

LWGT 

INTEGER*2 

5564 

M 

INTEGER*2 

33164 

MASK 

INTEGER*2 

5616 

N 

INTEGER*2 

33162 

NRNG 

INTEGER*2 

33 1 60 

OV 

REAL 

33090 

P 

REAL 

33492 

R 

REAL 

2 

R1 

REAL 

33026 

R2 

REAL 

33034 

RLPART 

real 

33503 

SI 

real 

1 03 1  a 

S2 

real 

21634 

SIN 

SORT 

TEMP 

real 

33276 

TEMPI 

real 

33298 

TRM 

real 

*32450 

U 

real 

33320 

VERTEX 

REAL 

11616 

XLOL 

REAL 

33133 

XR 

real 

33472 

XUPR 

real 

33146 

YLOL 

real 

33142 

YR 

real 

33476 

YUPR 

real 

33150 

I 

real 

0 

ZF 

real 

2704 

ZFMAX 

real 

33480 

ZFMIN 

real 

33434 

ZLEV 

real 

5408 

ZLOW 

real 

33154 

INTRINSIC 


M 1 crosoi 


-0?JTiwr,77 


/WORK  / 


/WORK  / 


/WORK  / 


INTRINSIC 

INTRINSIC 


/WORK  / 


/WORK  / 
/WORK  / 


/WORK  / 


V  ✓ 


□  Line#  1 


Microsoft  F0RTR0N77  V3.  EO  0£/3«- 


Name  Type 


MAIN 

ME3HS 

NUMBER 

P3D2D 

PLOT 

PLOTS 

SYMBOL 

WINDOW 

WORK 

ZCNTUR 

ZLEVEL 


PROGRAM 

SUBROUTINE 

SUBROUTINE 

SUBROUTINE 

SUBROUTINE 

SUBROUTINE 

SUBROUTINE 

SUBROUTINE 

COMMON 

SUBROUTINE 

SUBROUTINE 


Pass  Oni 


No  Errors  Detected 
448  Source  Lines 


•  V- V  VvVv 


V-  s’C •**',*' 


APPENDIX  E 


=*s» 

03-27-63 


3  Linas*  1  7 

1  *LARGE 

2  SSTORAGE :  3 

3  SBAGESirE:38 

4 


Microsof*  =CRTRAN77  73.30  03/34. 


5  C 

6  C 

7  C 
3  C 
3  C 

to  C 

11  c 

12  C 

13 

14  C 

15 

16 
17 
13 

13  C 
20 
21 
22 

23 

24 
23 
26 


* 

*  THE  PURPOSE  OF  THIS  PRGGRAM  IS  TO  CODE  THE  1-D  (DISCRETE 

*  TIME)  SYSTEM  TO  A  2-D  SPACIAL  SYSTEM. 

* 

*  EVANGELOS  THEOFILOU 

**************************************************************H 

PROGRAM  2D-DATA-FIELD 

*****  VARIABLE  DECLARATIONS  ***** 

REAL  R (23,  623)  ,  S(25,  623)  ,  R1 (2) , R2(2) , TRM(30, SO) , IV (SO) , 

*  IMGPORT 
CHARACTER* 1  ANSWER 

»**•  VARIABLE  DECLARATIONS  FOR  PL0T33  ***** 

CHARACTER*20  CTEXT 

COMMON  /WORK  /Z (26,  26)  ,  ZF(26,  26) , X (630) , Y (630) , ZLEV(26) , 

*  LDIG  <26) , LWGT (26) , MASK (3000) , VERTEX (16) 


DATA 


XLOL/O.  0/ ,  YLOL/O.  0/,  XUPR/8.  3/,  YUPR/7.  0/, 
ZLQW/ 1. 0E33/,  IPROJ/O/, NRNG/ 100/ 


27 

c 

main 

PROGRAM  *♦■ 

23 

23 

c 

*****  ASK  THE 

REQUIRED 

VALUES  FOR  THE  MODE! 

30 

31 

2  WRITE  (*,403) 

32 

READ  (*, *)  N 

33 

IF  <  (N  .  LT.  3 

)  .OR.  (N 

. GT. 23) )  GOTO  2 

34 

3  WRITE  (*,*04) 

■*?* 

READ  <*, *)  M 

36 

IF  ( (M  . LT.  2 

)  .OR.  (M 

.ST.  23))  GOTO  3 

37 

33 

WRITE  (*,401) 

33 

READ  (*,200) 

ANSWER 

40 

IF  ( (ANSWER 

.  EQ.  ’  Y’> 

.OR.  (ANSWER  .EO.  ’ 

41  c  ****  FILL  0’ s  THE  TWO  DIMENTIONAL  GRID  OF  CONTROL  POINTS  **** 

42  DO  36  I  »  1,26 

43  DO  36  Z  -  1,26 

44  Z(I,J)  »>'->.  O 

43  36  CONTINUE 

46 

47  C  ****  ENTER  VALUES  FOR  Y  MATRIX  **•*» 

46  DO  37  I  -  1,M*N 

49  WRITE  (*,408)  ’Y  (’,!,’>:  • 

30  READ  (*,  *)  R(l,  I) 

31  37  CONTINUE 


157 


'>3* 


7 

END  I" 

IF  ( (ANSWER  .  EQ 


Microsoft  FQRTRAN77  V3.  20  02/84 


36 

37 

38  C 

39 
80 
81 
62 
S3 

64  C 
63 
66 
87 

sa 

69 

70 

71 

7£  C 

73 

74 
73 

76 

77 
73 

79  C 

ao 

ai 

aa 

83 

84 

85 

86 
87 

aa 

89 

90 

91 
96 

93 

94 
93 

96 

97 

98 

99 
100 
101 
102  C 


IF  ((ANSWER  .EQ.  ’  Y’  >  .OR.  (ANSWER  . EQ.  •  y’  )  )  SOTO  4 

WRITE  (*,402) 

READ  (*,200)  ANSWER 

IF  ((ANSWER  . EQ.  ’ Y’  )  .OR.  (ANSWER  . EQ.  ’ y’  ) )  THEN 
•••*  INITIALIZE  THE  TRANSITION  MATRIX  **** 

DO  98  I  *  1,M+N 
DO  38  J  »  l,M-t-N 
TRM(I,J)  -0.0 

98  CONTINUE 

****  ENTER  VALUES  FOR  TRANSITION  MATRIX  **** 

DO  99  I  »  1,  M*N 
DO  99  J  *  1,M+N 

WRITE ( *,  407)  ’  T (’  ,  I, ’ , ’ , J, ’ )  s  ’ 

READ  (*,  *)  TRM ( I, J) 

99  CONTINUE 
END  IF 

*****  INITIALIZE  R  AND  S  ARRAYS  ***** 

DO  100  I  *  1,23 
DO  100  J  =  1,625 
R  <  I ,  J)  *0.0 
S<  I,  J)  »  0.  0 

100  continue 

*****  INITIALIZE  INPUT  vector  ***** 

DO  101  I  *  1,30 
IV  ( I )  *  0.  0 

101  CONTINUE 

WRITE  (*,211)  ’ ENTER  INITIAL  CONDITIONS  FOR  HORIZONTAL  R#’ 
DO  1 02  I  =  1 , M 

WRITE  <*, 405)  *  R’  ,  I,  ’  j  ’ 

READ  <*,  *)  R  ( I,  1 ) 

102  CONTINUE 

WRITE  (*,211)  ’ENTER  INITIAL  CONDITIONS  FOR  VERTICAL  S#’ 

DO  103  I  *  1 , N 

WRITE  (*, 405)  ’ S’ , I, ’ s  ’ 

READ  <*,  *)  S(I,  1) 

103  CONTINUE 

WRITE  (*,211)  ’ENTER  VALUES  FOR  THE  INPUT  VECTOR’ 

IV  (  1 )  »  1.0 

WRITE  (*,406)  ’ aO 1 :  ’ 

READ  (*, *)  IV (M*l ) 

IF  ((ANSWER  .EQ.  ’  Y» )  .OR.  (ANSWER  . EQ.  ’  y’  >  )  GOTO  5 
*****  INITIALIZE  TRANSITION  MATRIX  ***** 


'  V.yf  -Mii  <' 


D  Line# 
103 


Z„  Microsoft  FQRTRAN77  v; 

D0  104  I  =  1,£5 
DO  104  J  =  1 , 25 
TRM( I, J)  =  0.0 

104  CONTINUE 

WRITE  (*,£11)  'ENTER  ELEMENTS  OF  THE  TRANSITION  lA’RIX’ 
TRM(M+1,M+N)  =  — IV(M+1) 

WRITE  (*,406)  'alO:  > 

READ  (*, *)  TEMP 
TRM(1,M)  =  -TEMP 
TRM(l,n+N)  =  -1.0 
WRITE  (*,406)  'all:  ' 

READ  (*, *)  TEMP 

TRM(M+1,M)  =  TEMP  -  TRM ( 1 , M)  *  TRM < M+ 1 ,  M+N) 

DO  105  I  *  £, M 
TRM(I, I-l)  =1.0 

105  CONTINUE 


DO  106  I  =  2+M,  M+N 
TRM  (I,  I-l)  =  1.0 

106  CONTINUE 

3  U  =  1 .  0 

DO  107  I  =  1,N*M 
DO  108  J  =  1 , M+N 

IF  (J  .  LE.  M)  THEN 
DO  109  JJ  =  1,M+N 

IF  (JJ  .LE.  M)  R(J, I+l)  =  R(J.  i*i) 

*  R(JJ, I)*TRM(J, JJ) 

IF  (JJ  .GT.  M)  RCJ,t+l)  =  R(J,I*i)  ♦ 

.  *  S(JJ-M,  I) *TRM( J,  JJ) 

109  CONTINUE 

R  ( J,  I  1 )  =  R  <  J,  I  +  l )  *  I V  ( J )  *u 
END  IF 

IF  (J  .  GT.  M)  THEN 
DO  110  JJ  =  1,M+N 

IF  (JJ  . LE.  M)  S<J-M, I+l)  =  S(J— M, l+l)  + 

*  R ( JJ,  I ) *TRM ( J,  JJ) 
IF  (JJ  .  GT.  M)  S(J— M,  I  +  l)  =  S(J— M,  I  +  l)  + 

*  S ( JJ-M,  I ) *TRM ( J,  JJ) 

HO  CONTINUE 

S(J-M,  I  +  l)  =  S(J-M,  I  +  l)  +  I V  ( J )  *tj 
END  IF 

108  CONTINUE 
U  =  0.  0 

107  CONTINUE 

WRITE  (*,£ll)  >*****  INPUT  VECTOR  **++*• 

WRITE  (*,300)  ( I V ( I ) , I  =  1 , M+N ) 


o 


Lines*  1 
154 


'iicrosoft  ?aar^cN77  vs.  so  ot 
«R‘  ~  (*,  2iu  ******  -ransi-ion  ■matrix  ****** 
do  in  :  =  i , m*n 


1 

156 

WR I Tc  ( *, 300) 

1 

157 

WRITE  <  * ,  210) 

i 

159 

1  1  1 

CONTINUE 

153 

160 

WRITE  (*,211)  ’  ** 

161 

DO  112  I  =  1 ,  (*1*N 

l 

166 

WRITE  (*,300) 

1 

163 

112 

CONTINUE 

164 

165 

WRITE  (*,211)  ’ 

1 66 

DO  113  I  =  l , n*n 

* 

167 

WR I TE  ( *,  300 ) 

i 

169 

113 

CONTINUE 

163 

170  c 

171 

172 

173 

174 
173 

176 

177 
179 
173 

190 

191  C 


*•*•**  FILL  0*5  THE  TUO  DIMENTIONAL  3RID  OF  CONTROL  POINTS  **** 
DO  114  I  =  i,2S 
DO  114  J  =  1,26 
-0.0 

114  CONTINUE 

4  DO  115  I  «  l,M 

DO  115  J  »  1. N 

Z ( I, J)  -  R< 1, ( I-i ) *N+J) 

115  CONTINUE 


192 

DO  113  I  =  l. 

1 

193 

WRITE  (*, *) 

1 

134  U3 

CONTINUE 

195  C 

*****  OUTPUT 

196 

WRITE  (*,203) 

197 

WRITE  (*,212) 

199 

DO  116  I  =  l. 

1 

193 

WRITE  (*,30 

1 

130 

WRITE  (*,21 

1 

131  116 

CONTINUE 

132 

133 

WRITE  (*,213) 

134 

WRITE (*, 421 ) 

135 

READ  (*,200) 

136 

IF  (  (ANSWER 

137 

DO  117  I  »  l. 

1 

139 

X  <  I )  -  0.0 

1 

133 

Y  (  I )  »  0. 0 

1 

200  1 l 7 

201 

CONTINUE 

202 

DO  119  I  =»  1,1 

1 

203 

X  (  I)  *  I  » 

1 

204 

V  (  I )  -  R  (  1 , 

205 

113 

CONTINUE 

206 

207 

18 

WRITE  (*,415) 

20S 

READ  <*,  199)  CTEXT 

209 

WRITE  (*,451) 

210 

READ  (*,200)  ANSWER 

21 1 
212  C 

*****  INITIALIZE  PL0T33 

213 

IF  <  (ANSWER  .  ED.  ’  Y’  ) 

214 

CALL  PLOTS  (0,  0,  2) 

215 

ELSE 

216 

CALL  PLOTS (0,  99, 99) 

217 

END  IF 

SIS 

219 

CALL  PLOT  (1.0,  1.0,  -3) 

220 

CALL  SCALE ( X, S.O,M*N, l) 

22 1 

CALL  SCALE ( Y,  4.  0,  M*N,  1 ) 
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CALL  STAXIS(0. 20, 0.20,0.  111,0.  112,  1) 

CALL  AXIS  <0.  0,  0.  0,  ’  X  AXIS’  ,  -6,  o.  0,  0.  0,  X  (M*N*1  )  ,  X  (M*N+2>  ) 
CALL  AXIS <0.  0,  0. 0,  ’ Y  AX  IS’ ,  S,  4.  0, 90. O, Y (M*N+1 >  ,  Y (M*N+2)  ) 
CALL  LINE  <  X,  Y,  M*N,  1,0,0) 

CALL  PLOT<0.  0,0.  0,-3) 

CALL  SYMBOL (1.0,  S.  73,  0.  23,  CTEXT,  0.  0,  20 ) 

CALL  SYMBOL (6. 0, 6. 3, 0. 2,  ’  1-0  DATA  FIELD’ , 0.  0,  14) 


230 

c 

***** 

OUTPUT  THE  GRAPH  ***** 

231 

CALL 

PLOT  <0.  0,  0.  0,  999) 

232 

WRITE 

(*,  4iS) 

233 

READ 

(*,200)  ANSWER 

234 

if  ((Answer  ,sq.  ’y»i  .or. 

(Answer  . e 

2-iS 

236 

19 

WR I TE  <  *,  419) 

237 

READ 

(*,200)  ANSWER 

238 

239 

if  ((answer  .  ne.  ’yv>  .and. 

(ANSWER  . 

240 

C 

***** 

ASK  THE  PARAMETERS  FOR 

the  graph 

241 

20 

WRITE 

(*, 210) 

242 

WRITE 

(*, *)  ’ ****  ENTER 

PLOT 

243 

WRITE 

(*, 410) 

24.4 

READ 

(*, *)  AZIM 

245 

WRITE 

(*,4ll) 

246 

READ 

(*,*)  ELEV 

247 

WRITE 

(*, 413) 

246 

READ 

(*, *)  itrim 

249 

WRITS 

(*,  414) 

230 

READ 

(*,  *)  ID IV 

231 

WRITE 

(*,415) 

232 

READ 

(*, 199)  CTEXT 

£33 

WRITE 

(*,  451 ) 

£34 

READ 

(*,200)  ANSWER 

parameters 


y’  >  > 


ie#  1 

iZa  C 

:57 

>5a 

:s3 

■so 

LSI 

>63 

>63 

>64 

165  C 

>66 

>67 

>68 

>63  C 

>70 

z~r  i 

>73 

>73 

>74 

>75 

>76 

>77 

>78 

373 

380 

381  C 
383 

383 

384 

385 

386 


7  Micr'oso 

*•*■*•»■*  INITIALIZE  PL0T38  *•*■»■*♦ 

IF  ((ANSWER  .ED.  ’ v» )  .DR.  (ANSWER  .ED.  ’ 
CALL  PLOTS  (0,  0,  3) 

ELSE 

CALL  PLOTS (0,  33, 33) 

ENDIF 

CALL  WIMDQWIXLOL,  YL.CL,  XUPR,  YUPR) 


OR  ■  RAr-7^  V3.  3 
THEN 


*****  DRAW  THE  MESH  SURFACE  OF  THE  GRAPH  ***** 

CALL  MESHS  <Z,  36,  36,  N,  M,  AZIM,  ELEV,  0.  5,  0.  5,  3.  35,  6.  2,  IDIV,  0, 
3,  I  PRO J,  1, ZLOW, 3,  ITRIM,  MASK, VERTEX) 


*****  ANNOTATION  OF  THE  GRAPH  ***** 

CALL  SYMBOL  (1.0,  6.  75,  0.  35,  CTEXT,  0.  0,  £0) 

CALL  SYMBOL  (6.  O,  6.  5,  0.  3,  ’  3-D  DATA  FIELD’ ,  O.  O,  14) 
CALL  SYMBOL (5. 2, 0. 3, O. 3, ’ AZIMUTH:  ’,0.0,10) 

CALL  NUMBER  (333.  0,  333.  0,  0.  3,  AZIM,  0.  0,  3) 

CALL  SYMBOL  <2.  5,  0.  0,  0.  3,  ’  ELEVATION:  ’  ,  0.  0,  10) 

CALL  NUMBER  (333.  0,  333.  0,  O.  3,  ELEV,  0.  0,  3) 

DY  =  (Z(l,  D/30.  0)  *  ELEV 

CALL  P3D3D  (1.0,  1. 0,  Z  (  1,  1>-DY,  XR,  YR> 

CALL  SYMBOL  (  XR,  YR,  0.  35,  ’  *»  ,  0.  0,  1 ) 

CALL  SYMBOLd. 0,0.  1,0.3,  ’*  =  ORIGIN’ ,  0.  0,  10) 


*****  OUTPUT  THE  GRAPH  ***** 

CALL  PLOT  (0.  0,  0.  0,  393) 

WRITE  (*,416) 

READ  (*,300)  ANSWER 

IF  ((ANSWER  .  EQ.  ’  Y’  >  .OR.  (ANSWER  .ED.  ’y’>)  GOTO  30 


387  31  WRITE(*, 413) 

388  READ (*,300)  ANSWER 

333  IF  <  (ANSWER  .ED.  ’  Y*  )  .OR.  (ANSWER  .ED.  ’ y’  )  )  THEN 

390  C  ****  FILL  0’ s  THE  TWO  DIMENTIQNAL  GRID  OF  CONTROL  POINTS 

331  DO  133  I  =  1,36 

333  DO  133  J  =  1,36 

333  ZF(I,J>  =0.0 

334  133  CONTINUE 

335 

396  ZFMAX  =  -9. 3E30 

397  ZFMIN  =  3.  9E30 

338  DN  =  (N-D/3.  0 

339  DM  =  (M-D/3. 0 

300  P  =  3. 141533 

30 1  DO  1 33  MM  =  1 ,  M 

303  DO  133  NN  =  1,N 

303  RLPART  =  0. O 

304  IMGPART  =  O. O 

305  DO  134  L  =  1, M 

306  DO  134  K  =  1 , N 


**** 


OJ  OJ  OJ  OJ 


Lire# 

307 

30a 

309 

310 
31  1 

312 

313 

314 

315 
31S 

317 

318 

319 

320 

321 

322 


Rl(l)  = 
Rl(2)  = 
R2  ( 1  >  =■ 
R2 (£)  = 
RLPART 


I MG PART 


Microsoft  -CRTRAN77  03.  20 
COS ( — 2*P* ( L— 1 >  *  t MM— DM— 1 ) / M ) 

SIN  < -2*P* ( L— 1 ) * (MM-DM-l ) /M) 

COS  <- 2*P* <K-1 ) *(NN-DN-I > /N) 

SIN  <-2*P* (K-l > * (NN-DN-1 > /N) 

=  RLPART  -r  Z(L,  K)*(R1  (1)  *R2  <  1 ) 

— R 1  (2)  *R2  (2)  ) 

=>  IMGPART  +  z (L,  K> *(R1  ( 1 ) *R£ (2) 

+R1 (2) *R2 ( 1 ) ) 


345 

346 

347 

348 

349 

350  C 

351 


CONTINUE 
ZF (MM, NN) 


ZF(MM, NN)  =*  SORT ( RLPART  **2 
IF  (ZF(MM, NN)  .  GT.  ZFMAX) 
IF  <ZF(MM, NN)  .  LT.  ZFMIN) 
CONTINUE 


+  IMGPART**2) 
ZFMAX  =»  ZF  (MM,  NN) 
ZFMIN  *  ZF (MM, NN) 


*****  OUTPUT  THE  ZF  MATRIX  ***** 

WRITS  (*,205)  ’ ***  FOURIER  TRANSFORMATION  ’ , M, »  X  ’ , N, • 
WRITE  (*,212) 

DO  135  I  =*  1,1*1 

WRITE  (*,300)  (ZF(I,J>,  J  =  1,N> 

WRITE  (*,210) 

CONTINUE 
WRITE  (*,213) 


WRITE  <*,  419) 

READ  <*, 200)  ANSWER 

IF  ((ANSWER  . NE.  ’ Y’ )  .AND.  (ANSWER  . NE.  >y>))  3DT0  22 


WRITE 

WRITE 

WRITE 

READ 

WRITE 

READ 

WRITE 

READ 

WRITE 

READ 

WRITE 

READ 

WRITE 

READ 


ASK  THE 
(*, 210) 
(*,  *)  ’ 
<*, 410) 
(*,  *) 

(*, 411) 

<*,  *) 

(*,  413) 

(*,  *) 

(*, 414) 


PARAMETERS  FOR  THE  GRAPH 


***  ENTER  PLOT  PARAMETERS  **-*’ 


ANSWER 


*****  INITIALIZE  PLQTS8  ***** 

IF  ((ANSWER  .ED.  ’V’)  .OR.  (ANSWER  . EQ. 

CALL  PLOTS  <0,  0,  2) 

ELSE 

CALL  PLOTS <0,  99, 99) 

ENDIF 


WRITE  (*,4S0> 


*V  *.  %  *.  'A  AV.  LA 


V'  V  ■l" 


V-  * 


D  Lirie*  l 
353 
253 
260 
261 
362 
262 
264 
263 

1  266 
1  267 

1  263 

263 

370 

371 

372 

373  C 

374 
275 

376 

377  C 
373 
373 
230 
331 

382 
333 
234 

383 
286 
337 
383 
383 

330  C 

331 

332 

333 

334 

335  i 
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READ  (*,200)  ANSWER 

CALL  WINDOW (XLGL,  YLQL,  XUPR,  YUPR) 

IF  ((ANSWER  . EQ.  ’V’)  .OR.  (ANSWER  .  EQ.  ’  y’  > )  THEN 
DLEV  =  (ZFMAX-ZFMIN) /FLOAT <M) 

CALL  ZLEVEL ( ZF,  26,  26,  M,  N,  DLEV,  ZLEV,  N> 

DC  136  I  =  1 , N 
LDIG(I)  =*  2 
LWGT ( I )  =1 

CONTINUE 

CALL  ZCNTUR  (  ZF ,  26,  26,  M,  N,  0.  5,  0.  5,  3.  25,  6.  5,  ZLEV,  lDIG,  LWGT, 

N,  0.10,  10) 

CALL  SYMBOL (5. 5, 0. O,  0.  2,  ’  CONTOUR  MAP’  ,  0. 0,  1 1 ) 

ELSE 

*****  DRAW  THE  MESH  SURFACE  OF  THE  GRAPH  ***** 

CALL  MESHS (ZF,  26,  26,  M,  N,  AZIN,  ELEV, 0. 5, 0. 5, 3. 25, 6. 5,  ID  IV,  0, 

3,  I  PRO J ,  1,  ZLOW, 3,  ’TRIM, MASK, VERTEX) 

*****  ANNOTATION  OF  THE  GRAPH  ***** 

CALL  SYMBOL (5. 5, 0. 3, 0. 2,  ’  AZIMUTH:  ’,0.0,10) 

CALL  NUMBER (333. 0, 333.  O, 0.  2,  AZIM, 0.  0, 2) 

CALL  SYMBOL (5. 5, 0. 0,  0.  2,  ’  ELEVATION:  ’  ,  0.  0,  10) 

CALL  NUMBER  <333.  0,  333.  0,  0.  2, ELEV, 0. 0, 2) 

DY  =»  ( ZF  ( 1,  1 )  /30.  0)  *  ELEV 

CALL  P2D2D (1.0, 1. 0, ZF ( 1 , 1 ) -DY, XR, YR) 

CALL  SYMBOL (XR,  YR, 0.  25, ’  *’  , 0. 0,  1 ) 

CALL  SYMBOL ( 1. 0, 0.  1, O. 2, ’ *  =  ORIGIN’ , 0. 0,  10) 

END  IF 

CALL  SYMBOL ( 1. 0, 6.  75,  0. 23,  CTEXT, 0. 0, 20) 

CALL  SYMBOL  (6.  0,6.  3,0.  2,  ’  2-D  DFT'  ,0.0,7) 

*****  OUTPUT  THE  GRAPH  ***** 

CALL  PLOT (0. 0, 0.  0,  333) 

WRITE  (*,416) 

READ  (*,200)  ANSWER 

IF  ((ANSWER  .EQ.  ’ Y’  >  -OR.  (ANSWER  .  EQ.  ’  y’  >  )  GOTO  30 


336 

WRITE  (*,417) 

337 

READ  (*,200)  ANSWER 

338 

IF  (  (ANSWER  .EQ.  ’  Y’  )  .OR.  (ANSWER 

•  EQ.  '  y’ 

) )  GOTO 

r* 

333 

STOP 

400 

401 

193 

FORMAT (A20) 

402 

200 

FORMAT (A) 

403 

£05 

FORMAT  </,  13X,  A23,  12, A3,  12,  A3, /> 

404 

210 

FORMAT  < ) 

405 

£11 

FORMAT (/, SX , 60 A ) 

406 

£  1  £ 

FORMAT (/, 2X, ’ (AZIMUTH  320. 0) ’ , 46X, ’ 

(AZ IMUTH 

230. 0) ’ , 

/) 

407 

£12 

FORMAT  </, EX , ’  (AZIMUTH  050. 0) ’  ,  46X, ’ 

(AZIMUTH 

140. 0) ’ , 

/  ) 

408 

200 

FORMAT ( 10 (F7. 2, IX) ) 

iw! 


bv 


is 


IA\ 


$ 


ns#  1 

7 

409 

400 

FORMAT  <3X,  \> 

410 

40  1 

FORMAT!/, 5X, 

411 

402 

FORMAT (/, 3X, 

412 

403 

FORMAT (/, 5X, 

413 

404 

FORMAT!/, SX , 

414 

405 

FORMAT (5X,  Al, 

415 

406 

FORMAT !5X, A5, 

416 

407 

FORMAT !5X, A2, 

417 

408 

FORMAT !5X,  A2, 

413 

403 

FORMAT !3X,  AS 

413 

410 

FORMAT !/, 5X, 

420 

411 

FORMAT (/, 5X, 

421 

412 

FORMAT!/, 5X, 

422 

413 

FORMAT!/, 5X, 

4  23 

414 

-ORMAT ! / , 5X, 

424 

415 

FORMAT!/,  5X, 

425 

416 

FORMAT (/,  5X, 

426 

417 

FORMAT!/, 5X , 

427 

418 

FORMAT!/, 5X, 

428 

413 

FORMAT!/,  5X, 

423 

420 

FORMAT!/,  5X, 

430 

421 

FORMAT!/,  5X, 

431 

432 

451 

FORMAT (/,  5X, 
END 
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DO  YOU  WANT  TO  FILL  THE  Z  .’MATRIX  ?  ’,\) 

DO  YOU  WANT  TO  FILL  THE  TRANSITION  MATRIX 
COLUMNS  OF  OUTPUT  FRAME (N- It *25) :  ’,V> 
ROWS  OF  OUTPUT  FRAME ( M= 1 t o25 ) :  ’,\> 

I£,A£, \) 

\) 

12,  Al, 12, A3, \> 

13,  A3, \> 

\) 

AZIMUTH  <0,  O  to 


60.0  DEGREES): 
ELEVATION (30. O  to  -30.0  DEGREES):  ’,\) 
NUMBER  OF  SMOOTHINGS:  ’,\> 

TRIM <0=N0, l*Xs, 2=Ys) :  ’,\) 

2, A  OR  a  SUBGRIDS:  ’ , \) 

TITLE  OF  GRAPHtUP  TO  20  CHAR):  >,\) 

DO  YOU  WANT  TO  CHANGE  PARAMETERS’?  ’,\> 
DO  YOU  WANT  TO  REPEAT  THE  PROCESS"’ ' 

DO  YOU  WANT  FOURIER  TRANSFORMATION 
DO  YOU  WANT  TO  MAKE  GRAPH  ?  ’,\> 

DO  YOU  WANT  CONTOUR  MAP  ?  ’,\> 

DO  YOU  WANT  TO  DRAW  CARVE  ?  ’,\) 

SEND  GRAPH  TO  THE  PRINTER ( Y  or  N) : 


\) 


\) 


’  ,  \) 


Name 

Type 

Offset  P 

Class 

ANSWER 

CHAR* 1 

30 

AZIM 

REAL 

134 

CCS 

CTEXT 

CHAR*20 

174 

INTRINS 

IC 

DLEV 

REAL 

284 

DM 

REAL 

23/1 

DN 

real 

226 

DY 

REAL 

206 

ELEV 

REAL 

138 

FLOAT 

INTRINSIC 

I 

INTEGER*2 

32 

I D I V 

INTEGER*2 

2/14 

IMGPAR 

REAL 

258 

IPROJ 

INTEGE.R*2 

22 

ITRIM 

INTEGER*2 

202 

IV 

REAL 

0 

LARGE 

J 

I NTEGER*S 

34 

JJ 

INTEGER*2 

110 

K 

I NTEGER*2 

270 

L 

INTEGER*2 

262 

LDIG 

INTEGER*2 

10552 

/WORK 

/ 

LWGT 

INTEGER*^ 

10604 

/WORK 

/ 

M 

INTEGER*2 

28 

MASK 

INTEGER*2 

10656 

/WORK 

/ 

165 


05  U 


0  Line*  1  7 


MM 

INTEGER*2 

N 

INTEGER*2 

NN 

INTEGER*^ 

NRNG 

INTEGER*^ 

P 

REAL 

R 

REAL 

R1 

REAL 

R2 

REAL 

RLPGRT 

REAL 

S 

REAL 

SIN 

SORT 

TEMP 

REAL 

TRM 

REAL 

U 

REAL 

VERTEX 

REAL 

X 

REAL 

XLOL 

REAL 

XR 

REAL 

XUPR 

REAL 

Y 

REAL 

YLOL 

REAL 

YR 

REAL 

YUPR 

REAL 

Z 

REAL 

TC 

REAL 

ZFMGX 

REAL 

ZFMIN 

REAL 

ZLEV 

REAL 

ZLOW 

REAL 

Name  Type 

AXIS 

LINE 

MG  IN 

ME3HS 

NUMBER 

P3D2D 

PLOT 

PLOTS 

SCALE 

STPX  IS 

SYMBOL 

WINDOW 

WORK 

ZCNTUR 

ZLEVEL 


223 

26 

246 

24 

234 

0 

LARGE 

0 

LARGE 

a 

LARGE 

254 

0 

LARGE 

78 

0 

INTRINSIC 

INTRINSIC 

LARGE 

34 

16656 

/WORK  / 

5408 

/WORK  / 

2 

210 

10 

7328 

/WORK  / 

6 

214 

14 

0 

/WORK  / 

2704 

/WORK  / 

218 

1 044 a 

/WORK  / 

18 

Size 

Class 

SUBROUTINE 

SUBROUTINE 

PROGRAM 

SUBROUTINE 

SUBROUTINE 

SUBROUTINE 

SUBROUTINE 

16720 

subroutine 

SUBROUTINE 

SUBROUTINE 

subroutine 

subroutine 

COMMON 

SUBROUTINE 

subroutine 

!u 

Ay*'?*?*1 


D  Line#  1  7 

Pass  One  No  »*'r*ors  Detectea 

^32  Source  Lines 
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